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Abstract 


Low Reynolds number homogeneous turbulence undergoing low Mach 
number isotropic and one-dimensional compression has been simulated by 
numerically solving the Navier-Stokes equations. The numerical simula- 
tions were carried out on a CYBER 205 computer using a 64 x 64 x 64 
mesh. A spectral method was used for spatial differencing and the 
second-order Runge-Kutta method for time advancement. A variety of 
statistical information was extracted from the computed flow fields. 
These include three-dimensional energy and dissipation spectra, two- 
point velocity correlations, one-dimensional energy spectra, turbulent 
kinetic energy and its dissipation rate, integral length scales, Taylor 
microscales, and Kolmogorov length scale. 

It was found that the ratio of the turbulence time scale to the 
mean-flow time scale is an important parameter in these flows. When 
this ratio is large, the flow is Immediately affected by the mean strain 
in a manner similar to that predicted by rapid distortion theory. When 
this ratio is small, the flow retains the character of decaying iso- 
tropic turbulence Initially; only after the strain has been applied for 
a long period does the flow accumulate a significant reflection of the 
effect of mean strain. In these flows, the Kolmogorov length scale 
decreases rapidly with increasing total strain, due to the density 
increase that accompanies compression. 

Results from the simulated flow fields were used to test one-point- 
closure, two-equation turbulence models. The two-equation models per- 
form well only when the compression rate is small compared to the eddy 
turn-over rate. A new one-point-closure, three-equation turbulence 
model which accounts for the effect of compression is proposed. The new 
model accurately calculates four types of flows (Isotropic decay, iso- 
tropic compression, one-dimensional compression, and axisymmetric expan- 
sion flows) for a wide range of strain rates. 
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Chapter I 
INTRODUCTION 


Turbulent flow undergoing compression is of great technological 
Interest; understanding it is therefore of considerable importance. For 
example, the compression stroke in an internal combustion engine mod- 
ifies the turbulence in the cylinder in important ways and thereby 
influences the nature of the combustion process and overall engine 
performance. Compression of turbulence is also important in shock wave- 
turbulent boundary layer interaction and other applications. However, 
current turbulence models, Including the popular k-e model, do not 
predict compression effects accurately. This is evidenced by their poor 
performance in piston-engine flow calculations. Reynolds (1980) and 
Morel and Mansour (1982) pointed out defects in the existing model and 
proposed new model constants. Although the predictions of the modified 
model look reasonable, they have not been tested against experimental 
data, thus leaving the model user in an uncomfortable position. 

1.1 Compression Effects on Turbulence in Engines 

The effects of compression on turbulence have been investigated 
both theoretically and experimentally. Rapid distortion theory (RDT) 
predicts the evolution of turbulence undergoing extremely rapid com- 
pression. Under the assumption that the strain is applied rapidly, 
nonlinear turbulence interactions are unimportant and the dissipation 
process is too slow to be significant. Thus, RDT is a linear, inviscid 
analysis of turbulence. It was originated by Taylor (1935) and was 
extensively developed by Batchelor and Proudman (1954). 

Hoult and Wong (1980) applied RDT to the prediction of the evolu- 
tion of the turbulence undergoing rapid compressions Computations based 
on this theory were compared with the available hot-wire measurements 
(Witze, 1977; and Lancaster, 1976), and reasonably good agreement was 
obtained . 

Hunt (1978) reviewed RDT and applied it to various types of flows. 
He found that in rapid, one-dimensional compression the vorticity that 
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is aligned with the compression axis does not change, while the vortic- 
ity in the other two directions increases. He also worked out the 
evolution of the components of the turbulence kinetic energy for the 
rapid-compression case. RDT describes a limiting case; the conditions 
for its validity are not often met in practice. However, it does pro- 
vide a useful analytical solution against which models and simulations 
can be tested. 

Semenov (1958) measured turbulence in a motored, single-cylinder 
engine with a hot-wire anemometer. He found that the turbulence 
intensity at top dead center (TDC) was affected primarily by engine 
speed. Winsor and Patterson (1973) used a hot-wire anemometer to 
measure velocity fluctuations in a motored engine and found that the 
turbulence decayed during compression. Lancaster (1976) used a tri- 
axial, hot-wire anemometer to study turbulence in a motored single- 
cylinder engine; the measured turbulence intensity was found to decrease 
quite rapidly, reaching a minimum near TDC. He attributed this result 
to changes in flow orientation or pattern and viscous dissipation. He 
also found that the spectrum of the turbulence energy shifted to lower 
frequencies during compression. These results supported his argument 
that turbulence decays during compression. The hot-wire experiments of 
Witze (1977) were performed on a motored, single-cylinder engine. He 
found that turbulence intensity increases during the compression stroke, 
but later (1980) admitted that his data-reductlon techniques overestima- 
ted the turbulence intensity. The problem stems from the fact that a 
hot wire does not measure the velocity directly, but responds to mass 
flux. As a consequence, it is necessary to account for the thermody- 
namic state of the gas when reducing the anemometer signal. There is 
reason to suspect that this problem is present in much of the published 
data taken with hot wires in engines. 

More recently, investigators have used laser-Doppler anemometers 
(LDAs) to study in-cylinder fluid mechanics in internal-combustion 
engines. The group at Imperial College did a series of investigations 
on this subject, Morse, Whitelaw, and Yianneskis (1979, 1980) made LDA 
turbulence measurements in motored, piston-cylinder assemblies without 
compression; only the intake and exhaust strokes were allowed in their 
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apparatus. Arcoumanis, Bicen, and Whitelaw (1982) studied turbulence in 
a model four-stroke motored engine with a compression ratio of 3.5. 
They concluded that the turbulence level decays considerably in the 
early part of the compression stroke, while in the second half of the 
stroke the piston configuration becomes increasingly Important. A flat 
piston produces an almost constant turbulence level. A cylindrical bowl 
has no significant effect on the flow, independent of the bowl depth 
within the practical range. However, the addition of a lip causes a 
substantial inward motion of the air toward the bowl but no apparent 
effect on the turbulence levels outside the bowl. The same authors 
(1983) also studied the squish effect on turbulence near TOC of com- 
pression in a motored model internal-combustion engine with compression 
ratio 6.7. They found that the squish did not alter the overall turbu- 
lence levels, which were comparable to those obtained with the flat 
piston. They claimed that squish is an important contributor to the 
mean motion rather than to the turbulence. 

At Stanford, Richman (1982) developed a Flow Diagnostics Engine 
(FDE) which is a single-cylinder engine with a transparent cylinder made 
from single-crystal sapphire. The valve motion can be completely con- 
trolled by a minicomputer. A special Schlleren system was developed to 
visualize flows within the engine cylinder. His visualizations sugges- 
ted that the scale of the turbulent motions decreases as the flow is 
compressed . 

To model the effect of compression on turbulence, Watkins (1977) 
extended the k-e model for incompressible flows developed by Launder 
and Spalding (1974), by adding a dilatation term to the modeled dissipa- 
tion equation. The new term had a constant equal to unity, a natural 
outcome of the derivation. Gosman and Watkins (1977), Gosman and Johns 
(1978, 1980), Gosman, Johns, and Watkins (1980) all employed this 
extended k-e turbulence model to predict the flow in reciprocating 
engines. Grasso and Bracco (1983) also introduced this k-e model into 
a modified version of a code developed at Los Alamos by Butler et al. 
(1979), and used it to study the squish effect. Ahmadi-Bef rui , Gosman, 
Lockwood, and Watkins (1981) made calculations with three different 
values of the model constant: 0, 1.0, and -0.373. Gosman and Harvey 
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(1982) used these constants to predict the behavior of a Diesel engine. 
Gosman and Jahanbakhsh (1981) reported that the predictions of turbu- 
lence intensity using -0.373 as the model constant are marginally better 
than those found using the alternative values. Ramos, Humphrey, and 
Sirignano (1979) used a purely "ad-hoc" turbulence model. 

Reynolds (1980) pointed out that Watkins* model predicts growth of 
the length scale during compression independent of the rate of compres- 
sion. Morel and Mansour (1982) found that, in an engine simulation, the 
length scale predicted by the Watkins model reached values several times 
the cylinder height near TOC; this is clearly incorrect. Reynolds 
(1980) applied rapid-distortion analysis to the isotropic compression 
case and suggested that the model constant in the dilatation term be 
-2/3. Borgnakke, Arpaci, and Tabaczynski (1980) independently found the 
same value. Morel and Mansour (1982) extended Reynolds* work to encom- 
pass three different "modes" of compression, i.e., isotropic, one- 
dimensional, and two-dimensional compressions. Ibis led to a different 
model constant for each case. Mansour (1982) employed this modified 
k-e model to study dispersion in stratified-charge engine cylinders 
during the compression stroke. His results are model-dependent and, in 
the absence of accurate data, they cannot be relied upon. 

El Tahry (1983) derived exact transport equations for turbulence 
kinetic energy and its dissipation rate. An order-of-magnitude analysis 
was applied to the terms in these equations, and the small terms were 
neglected. It turned out that the k-equation retains its usual form, 
while four new terms are added to the dissipation equation; the constant 
in the dilatation term was found to be -1/3. Results of calculations 
using this modified dissipation equation are more physically plausible 
than those obtained with Watkins' version of the model. However, this 
does not answer the question of whether the present version of the 
k-e model is adequate for use in engine flow computations . 

Dussauge, Gaviglio, and Favre (1978) suggested that production of 
turbulence can be broken into two parts: (i) a "dilatation" production 
and (ii) an "isovolumetric” production related to the incompressible 
component of the strain. This suggestion will be used in the present 
work. 
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In summary, RDT provides a limiting solution demonstrating how com- 
pression affects turbulence. In attempts to understand the effects of 
compression on turbulence for the range of parameters encountered in 
engineering flows, people have looked to experiments. Due to the 
complexity of the flows studied and the difficulties associated with 
measurements in complex flows, only qualitative results have been ob- 
tained. Model builders have suggested various turbulence models aimed 
at representing the effects of compression; these models differ, espe- 
cially with respect to the constant associated with the dilatation 
term. Users have employed these models in codes aimed at predicting 
fluid motion Inside internal-combustion engines with only moderate 
success. It is difficult to identify the source of the discrepancy 
between the numerical predictions and the experimental results. A com- 
mon conclusion by model-users is that the turbulence models appear to 
predict the general behavior of the turbulence field. However, the 
quality of experimental data does not allow accurate evaluation of the 
models . 

1 .2 Turbulent Flow Simulation 

Newtonian fluid motion is governed by the Navier-Stokes equations . 
However, the wide range of length scales present in turbulent flows 
makes simulation of them via solution of the Navier-Stokes equations 
impossible at the Reynolds numbers of interest in engineering flows. 
However, full numerical simulation of turbulent flows at low Reynolds 
number in simple geometries has recently become feasible. In this 
approach, the only errors made are numerical ones which can be con- 
trolled. These simulations can be regarded as numerical experiments, 
and the results complement laboratory measurements. Such simulations 
can provide understanding of turbulence phenomena. A unique advantage 
of numerical simulation is that, with it, one can single out a particu- 
lar effect and study its influence on turbulence. Full simulations are 
currently limited to simple, low-Reynolds-number turbulent flows and 
require supercomputers; they are impractical for direct engineering 
applications. However, they can be used as a research tool for inves- 
tigating both the physics of turbulence and the models used to represent 
it. 
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We limit our discussion of full simulation to homogeneous flows. 
In these flows, the statistical turbulence quantities are independent of 
spatial position; they develop in time. In homogeneous flows without 
mean strain or shear, the turbulence decays with time; when mean strain 
or shear are applied, the kinetic energy of the turbulence may increase 
with time. 

Homogeneous isotropic turbulence has been simulated by Orszag and 
Patterson (1972), Schumann and Patterson (1978), Clark, Ferziger, and 
Reynolds (1979), and many other people. The results have been used to 
determine turbulence model constants associated with dissipation, and it 
is usually the first flow simulated by people doing full simulations. 

The next group of flows contains those in which there is energy 
exchange between the various components of the Reynolds stress (redis- 
tribution) in addition to dissipation, but no direct production of 
turbulence energy . There are two such flows : homogeneous turbulence 
with rotation (which has been simulated by Rogallo , 1981, and Bardina, 
Ferziger, and Rogallo, 1985) and homogeneous flow relaxing from the 
effects of strain or shear. Schumann and Herring (1976), Schumann and 
Patterson (1978), and Rogallo (1981) simulated turbulence undergoing 
relaxation from axisymmetric contraction strain. Lee and Reynolds 
(1985) studied the relaxation of turbulence after various types of 
strain; some details of their simulations will be given later. 

Hie final group of homogeneous flows contains flows in which all of 
the phenomena that are possible in homogeneous flows actually occur: 
production, dissipation, and redistribution. There include homogeneous 
sheared turbulence and various types of homogeneous strained turbulence. 

Full simulation of homogeneous sheared turbulence was performed by 
Rogallo (1981). Shirani, Ferziger, and Reynolds (1981) added a passive 
scalar to homogeneous shear flow, and Feiereisen, Reynolds'; and Ferziger 
(1981) studied the effects of compressibility on homogeneous sheared 
turbulence . 

Full simulation of homogeneous strained turbulence was performed by 
Rogallo (1981). He considered the axisymmetric contraction and plane 
strain cases. Dang (1985) simulated homogeneous turbulence subjected to 
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two successive plane strains. Lee and Reynolds (1985) did the most com- 
prehensive study of homogeneous turbulence subjected to plane, axisym- 
metric contraction, axisymmetric expansion strains, and combinations of 
them. 


1.3 Motivation and Objectives 

Turbulence undergoing compression is of considerable technological 
Interest, but the effects of compression on turbulence are not well 
understood. RDT provides solutions for cases which are far from those 
that arise in engineering applications. Due to the limitations of 
measurement techniques, experiments have provided only a qualitative 
picture of the effects of compression. As the state-of-art of exper- 
imental techniques advances, it will shed more light on this subject. 
On the other hand , model-builders have developed various turbulence 
models intended to represent the effects of compression but have insuf- 
ficient data with which to validate the models. 

Full simulation of turbulent flows is a research tool which allows 
one to better understand the physics of turbulence and is a valuable 
supplement to laboratory measurements. In particular, full simulation 
of turbulence undergoing compression is a unique tool that can be ap- 
plied to answer some of the questions raised above. In this work, we 
perform full simulations of homogeneous turbulence undergoing isotropic 
and one-dimensional compression and use the resulting flow fields as a 
data base for turbulence model development. The simulations were per- 
formed on a CYBER 205 computer. 

Our objectives are : 

1 . To simulate homogeneous turbulence subj ected to isotropic and one- 
dimensional compression at low Reynolds number by solving the 
Navier-Stokes equations numerically. 

2. To study the effects of compression on turbulence using the data 
generated by full simulations . 

3. To test various turbulence models for compressed flow by comparing 
them with results of the full simulations . 
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4. To construct a turbulence model which accurately predicts compres- 
sion effects on turbulence. 

In Chapter II we present the mathematical foundations of these sim- 
ulations. The numerical methods employed are discussed in Chapter III. 
In Chapter IV we describe the data management and the performance of the 
code on a CYBER 205 computer. In Chapters V and VI we present the re- 
sults of isotropic and one-dimensional compression simulations, respec- 
tively. Turbulence model testing and a new model are presented in 
Chapter VII. Chapter VIII contains the conclusions. 
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Chapter II 

MATHEMATICAL FORMULATION 


This chapter describes the flows simulated and presents the mathe- 
matical formulation including the governing system of equations, the 
solution for the inviscid flow, and the coordinate transformation used 
to transform the governing equations from a Cartesian coordinate system 
to the Lagrangian coordinate system in which the simulations were per- 
formed. 

2.1 Problem Statement 

In this section we describe the flows to be simulated numerically. 
Since the principal objective is to understand the effects of compres- 
sion on turbulence, we choose a flow to study which is not affected by 
the presence of solid boundaries. Homogeneous turbulence offers this 
possibility. In particular, we shall consider the case of a fixed mass 
of turbulent fluid contained within a rectangular parallelepiped, the 
opposing sides of which can move inward or outward with time. The tur- 
bulence may undergo uniform isotropic compression if all three pairs of 
sides move Inward at same rate; this case is illustrated in Figure 2-1. 
This case simulates the squish effect of an engine which has a cup-in- 
piston design. If only one pair of sides moves inward as shown in Fig- 
ure 2-2, the turbulence undergoes one-dimensional compression. This 
case simulates the compression stroke in an internal combustion engine 
with a flat piston and is also related to the compression of turbulence 
passing through a shock wave. 

There are several Important assumptions on which the simulations 
are based. These are stated below. 

• The fluid is Newtonian and obeys the ideal gas law. 

• Body forces are negligible. 

• The compression is adiabatic. 

• The Mach number is sufficiently small that sound waves play no 
significant role. Under this assumption, the fluid is com- 
pressed so that the fluid density depends only on time, not on 
space . 
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• Since there is no heat release and the Mach number is small, the 
fluctuations of temperature can be neglected which means that we 
can assume the fluid properties to be functions of time alone. 


2.2 Mathematical Formulation 

The basic equations of motion for a Newtonian fluid are the Navier- 
Stokes equations : 


M 


5t + 0> u i),i “ 0 

»(p u i) . , , 

3t + (e u i u j),j " a ij,j 


’U 


2 

— Pfi + ufu + U — — U A 1 

P ij 1,3 l,i 3 k,k °ij j 


(2-1) 

(2-2) 

(2-3) 


where t, p, u, o» p, 


and 


are time, density, velocity, stress, 


pressure, and absolute viscosity, respectively. Indices run from 1 to 
3, and a repeated index in any term implies summation. Subscripts after 
commas denote partial differentiation, e.g., u^ j ■ du^/dxj . 

All flow quantities may be decomposed into mean and turbulence 
components : 


u^x.t) - U (x,t) + u^(x,t) 


(2-4) 


p(x,t) « p(t) + p*(x,t) 


(2-5) 


p(t) - p(t) 


(2-6) 


Substituting these into the equations of motion, (2-1), (2-2), and 

(2-3), taking the ensemble average and applying the assumption of homo- 
geneity of the turbulence so that (u'u’ ) » 0, we obtain the following 

1 3 ,3 

equations for the mean components: 


f? + t’V,i * 0 


(2-7) 


3(P V 

- 0 


(2-8) 
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Subtracting the mean equations from the total equations, we obtain the 
basic equations of motion for the turbulent component of the flow field. 



*(p u J) 

3t 



+ y; + 


(2-9) 



V U 




( 2 - 10 ) 


2.3 Mean Flow 

The governing equations of the mean flow are Eqs. (2-7) and (2-8). 
The mean flow provides the driving force for the turbulence. Eq. (2-8) 
shows the mean flow is driven by its inertial force because the right- 
and side of the equation which represents the surface forces effect is 
zero. 

Due to the assumption of homogeneity of the turbulence, the mean 
velocity field must be linear in the spatial coordinates 


U i,/j 


( 2 - 11 ) 


where U 
has the form 


depends only on time. 


The mean-deformation tensor 




(2-12) 


In isotropic compression, Sj(t) - S^Xt) ■ 83 (c); S£(t) * 83 (c) ■ 0 in 
the one-dimensional compression case. S i (t) is V /X p (t), where V p 
is the compression speed (assumed to be constant in all runs herein) and 
X p (t) is the box length at time t. The mean velocity vanishes at one 
end of the box while at the other end of the box it is equal to the com- 
pression speed V . The mean velocity profile for one-dimensional com- 
pression case is illustrated in Figure 2-3. 


The density of the fluid can be obtained by applying the mass con- 
servation law. For isotropic compression cases 
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(2-13) 


p(t) Xp(t) - p(0) 1 ? 

while, for the one-dimensional compression cases: 


p(t) x (t) ” p(0) L q 


(2-14) 


where L Q ,j enote8 the Initial box length. 

Due to the assumptions of adiabatic compression of the mean flow 
and the ratio of specfic heats, y, of the fluid being constant, the 
mean temperature can be obtained from thermodynamics : 


T(t) 

T(0) 



(2-15) 


The mean pressure can be obtained from the assumption of ideal gas beha- 
vior of the fluid, i.e.. 


P(t) - p(t) RT(t) 


(2-16) 


where R is the gas constant. 

The absolute viscosity of the fluid y(T) is assumed to have the 
following functional form in the range of the mean temperature we con- 
sider (Toulouklan et al., 1975): 


y(T) 


145.8 T 1 * 5 _-8 

x 10 


kg/(m-sec) 


110.4 + T 


for 250K < T < 850K (2-17) 


Since the mean temperature is a function of time only, we can rewrite 
Eq. (2-17) as 


u(t) m /T(t)\ 1<5 / 110.4 + K0) \ (2-18) 

V (0) \ T(0) / \ 110 .4 + T( t ) / 

Hie second term on the right hand side of Eq. (2-18) can be approximated 
by: 


! 110.4 + T(0) \ / T(t) 

\ 110.4 + *T( t ) / \ T(0) 


(2-19) 
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Substituting Eq.(2-19) into (2-18), we get the absolute viscosity as a 
function of time 


u(t) 

W<0) 


. ,0.75 


/ T(t) V 
\ T(0) / 


(; 


We summarize the mean flow for each type of compression as follows : 
Isotropic compression; 


UjCx.t) - 

V 

P x 

X (t) X 1 

p 

u <x.t) - 

V 

p x 

X (t) X 2 
P 

U 3 (x,t) - 

V 

p x 

X (t) 3 

P 

7(0 - 7(o) 


T(t) - T(0) 

/ L vl.2 

( ° \ 

lyW 

p(t) ■ 7(0) 



y(t) - y(0) 

( L o \ 0 - 9 

{ V‘>; 


C 


:- 20 ) 


:-2i) 
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One-dimensional compression : 


UjCx.t) 


v 

P_ 

x (t) 
P 


*1 


U 2 (x,t) - 0 


U 3 (x,t) ■* 0 

__ __ L 

p(t) * p(0) (2-22) 

P 



These mean properties will be used in the simulations. 

2.4 Coordinate Transformations 

The equations of motion for the turbulence field in the inertial 
laboratory frame are Eqs. (2—9) and (2—10). Periodic boundary conditions 
will be Imposed in this coordinate system. Although these boundary con- 
ditions are not exact, they are the best that one can do in simulating 
homogeneous flows and allow the turbulence to develop with a minimum of 
external interference. Previous simulations (Rogallo, 1981; Feiereisen 
et al., 1981; Shirani et al., 1981) showed the usefulness of periodic 
boundary conditions for simulating homogeneous turbulence. 

The presence of in the convective term in Eq. (2-10) makes 

this term aperiodic in space and prohibits the application of periodic 


14 



boundary conditions to the equation as a whole. To overcome this prob- 
lem, we need to introduce a coordinate transformation which eliminates 
the non-periodic terms. This new spatial coordinate system moves with 
the mean flow; it is essentially a Lagrangian coordinate system i.e. we 
are riding on the mean flow watching the evolution of turbulence. 

We define this transformation for the general case by 


Vj 


t r 


(2-23) 


where is a tensor relating the transformed coordinates to the 
Cartesian coordinates, x^ and t' are the transformed coordinates, and 
x^ and t are the Cartesian coordinates. Following Rogallo (1981), 
the transformation tensor Bjj is the solution of the following set of 
ordinary differential equations : 



dt 


B V 
ik k,j 


(2-24) 


subject to the initial conditions: 


ha ' s y • 

For the isotropic compression case, 
is 


at t ■ 0 

the mean-deformation tensor 


(2-25) 



and 


where Vp 



Bj j(t) * exp 
is a constant and 


V /X 0 

p p 

0 V x „ 

p p 

0 0 

H * 

Xp is equal to 



V° 


h> * V P *• 


2^) ™ B33 ( t ) * Bjj(t) 


(2-26) 


(2-27) 


(2-28) 
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Eq.(2-23) then becomes 


x f 

X 1 




ytT x 2 


ITTtT * 3 


(2-29) 


For the one-dimensional compression case, the mean-deformation tensor 


d i,j 16 


U 


i.j 


so 


V /x 

p p 

0 

0 


0 

0 

0 



(2-30) 


*u w ■ 




P 

B33 ( t ) 


o 

x_TtT 


(2-31) 

(2-32) 


Equation (2-23) then becomes 


t! - 


jTTtT 1 

P 


(2-33) 
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Figure 2-4 shows how the new coordinate system moves with the mean flow. 
Transformation of the equations from the Cartesian coordinate system to 
the new coordinate system is done by using the chain rule. The result- 
ing equations will be given below. 


2.5 Isotropically Compressed Turbulence 

The governing equations for the isotropic compression cases in the 
Lagrangian coordinate system are obtained using Eq. (2-29) and chain 
rule . We get 


* v x: 

9 _ _ pl 

dt x p (t) ax 


a 

L o a 

3Xj 

x <t) axT 
p i 

a 

L 

o a 

ax 2 

x (t) ax' 
P 2 

a 

L o a 

ax 3 

x ct) axT 
p 3 

a _ 

v x' 

p 2 a _ 


v x: 


x (t) ax' 


TtT axi + Tt r 


(2-34) 


a 2 . 

* 

3 2 

ax^Xj 

Xp(t) 

ax|axj 

3 2 . 

L 2 

o 

3 2 

ax 2 ax 2 

x 2 (t) 

ax’ax^ 

3 2 . 

L 2 

o 

a 2 

ax 3 ax 3 

x p(t) 

ax^ax^ 


Substituting Eqs. (2-34) and (2-21) into (2-9) and (2-10), we get the 
transformed continuity equation 


auj 

ax^ 


(2-35) 
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The momentum equations become 



9t 


+ 



+ 






a *j a *i 


(2-36) 


where Vp, L Q , and Xp are the compression speed, initial box length, 
and the box length at time t respectively. 


Density and absolute viscosity can be related to the instantaneous 
box length by: 



Since the troublesome convection terms has been transformed away, the 
transformed equations admit periodic boundary conditions. 


2.6 One-Dimens ionally Compressed Turbulence 

The governing equations for one-dimensional compression cases in 
the Lagranglan coordinate system are obtained in much the same way as 
for the isotropic compression case. However, because the mean flow in 
these two cases are different, the transformed governing equations are 
not the same. We get: 

L 

3 . o 3 

3 x x (t) ixT 

1 p 1 


3 _ _3_ 

3x 2 3x£ 



3t 


v x] 

— E-J: — § + 

x p 3XJ 



(2-37) 
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9X38X3 



3X33*3 


Substituting Eqs.(2-37) and (2-22) into (2-9) and (2-10), we get the 
transformed continuity equation 



(2-38) 



where Vp, L Q , and x p are the compression speed, initial box length, 
and the box length at time £ respectively. 


Density and absolute viscosity can be related to the instantaneous 
box length by: 

p(t) - 7(0) 

/L \0.3 

y(t) - y(0)^| 

There are no known analytical solutions of the equations of motion for 
the turbulence field, Eqs. (2-35), (2-36), (2-38), and (2-39). Numeri- 
cal approximations are needed to solve these equations. The simulations 
will be carried out in a finite domain with specified boundary and ini- 
tial conditions. These are described in the next chapter. 
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Fig. 2-2. Schematic description of homogeneous, 
one-dimensional compression flow. 


20 



X 2 

4 



Fig. 2-3. The mean velocity profile for one- 
dimensional compression flow. 




pig. 2-4. Lagrangian coordinate system for one- 
dimensional compression simulations. 
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Chapter III 
NUMERICAL METHOD 

In this chapter full simulation of turbulent flow is introduced, 
numerical methods used to perform the simulations are presented and 
techniques for removing aliasing error and method to generate initial 
velocity field are described. 

3.1 Introduction 

Since the equations of motion for the turbulence are highly non- 
linear, they are not amenable to classical analytical approaches. 
Numerical approximations are needed to solve these equations. Detailed 
simulations of turbulent flows can be used to both generate physical 
understanding and to improve the models. Full-turbulence simulation is 
the numerical solution of the exact Navler-Stokes equations. The only 
errors made are numerical ones which can be kept within a desired toler- 
ance. 

In full- turbulence simulations, the number of spatial grid points 
is determined by two constraints : the computational box has to be bigger 
than the size of the large eddies if the simulation is to capture them. 
On the other hand, the mesh size has to be smaller than the smallest 
eddies in order to resolve them. 

Ihe smallest scale of Importance is the Kolmogorov length scale 
which is defined as : 


n " [v 3 / (3-1) 

where v is the kinematic viscosity and e is the dissipation rate. 
Physically q is the size of the eddies for which viscosity is impor- 
tant. 

The scale of the energy-containing eddies is usually defined by 

l » q 3 /e (3-2) 
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where q is turbulent velocity. Therefore, the number of mesh points 
required to resolve both the large and small scales along one direction 
has to be roughly 

A - (St)'* - Ref 4 (3-3) 

At a Reynolds number based on turbulence .quantities of 10,000, a 
mesh 1000 points on a side is required. Clearly, this is beyond the 
capacity of any existing computer. Full-turbulence simulations are 
currently limited by computer resource availability to low turbulent 
Reynolds numbers. In our case, with 64 x 64 x 64 mesh points, the 
Reynolds number based on q and I has to be less than about 250. 

3.2 Approximation of Spatial Derivatives 

Spectral methods are often used in turbulent flow simulations, 
because, for sufficiently smooth fields, they are very accurate. In 
addition, the periodic boundary conditions used in homogeneous turbu- 
lence simulation make spectral methods based on Fourier expansions quite 
natural and easy to apply. We shall use the pseudo-spectral method to 
compute the spatial derivatives in terms of the data at mesh points. 
Use of this method in the simulation of turbulence was pioneered by 
Orszag and Patterson (1972) and used by many others since then. This 
method is made efficient by the Fast Fourier Transform (FFT) algorithm 
developed by Cooley and Tukey (1965) which is most efficient when the 
number of mesh points in each direction is a power of 2. Since this 
method treats each direction independently, we shall consider only the 
one-dimensional case here. 

A. Spectral Method 

To see how the method works, consider any function F(x) periodic 
with period L. The values of such functions at equally spaced mesh 
points can be represented in terms of a discrete Fourier series as 
follows s 

B/2 ' 1 - Ik X, 

F(x.) - £ F(k n > e n 3 (3-4) 

n»-N/2 
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where N is the number of mesh points and 
xj - hj, j - 0, 1, 2, 3, ..... N-l 

h ■ L/N is the mesh size 


2un/L, n 


-N/2, 


,, 0, N/2 - 1 


The coefficients of the Fourier series can be obtained from the inverse 
of Eq. (3-4) : 


N-l 


F(k ) - ~ IT F(x ^ e 

n N f— ' j 

j=o 


-ik x 

n j 


(3-5) 




Equation (3-4) can be rendered an interpolation scheme by replacing 
with a continuous variable x. Then the derivatives of F(x) with 


respect to x at 
polant : 


Xj can be approximated by differentiating the inter- 


dF(x ) 

N/2-1 

ik x. 

j 

■ E 

n»-N/ 2 

A n 1 

ik F(k ) e J 

n n 

dx 

d 2 F(x ) 

N/2-1 

• ik x 

j 

■ E, 

-k 2 F(k ) e n i 
n n 

2 

dx 

n— N/2 



(3-6) 


(3-7) 


The results are an extremely accurate estimates of the derivatives. 

The nonlinear terms in the equations of motion (2-36) and (2-39) 
introduce the possibility of aliasing errors into the numerical simu- 
lations. To understand the problem, let us consider the one-dimensional 
case with period L ■ 2* - Nh. If two functions represented as discrete 
Fourier series of the form (3-4) are multiplied, we obtain: 

N/2-1 N/2-1 . N-2 




a(Xj)b(x^) 


2 L e 

n*-N/2 m*-N/2 


i(n+m)x. 


2 M 

8—N 


isx. 


N/2-1 isx -N/2-1+N 


isx. 


N-2-N 


isx. 


2 «.). 2 «s). j + 2 *•>• 3 

s— N/2 s=-N+N s-N/2-N 


Aliasing Error 


(3-8) 
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Hie product series contains wave-numbers in the range -N < s < N-2 , 
which includes wavenumbers outside the range supported by the grid we 
are using. Fourier functions with wave-numbers outside the domain 
-N/2 jC s N/2 - 1 are interpreted as waves that belong to the computa- 
tional domain. On the grid, the functions are identical, but their 
derivatives are not identical. This introduces an error called alias- 
ing. The alias error resulting from nonlinear terms can be removed by 
two different methods (Patterson and Orszag, 1971) which will be ex- 
plained below. 


B. Alias Removal By Truncation 

One alias-removal method is the "2/3 rule". It consists of the 
following : 

(1) Given a(n) and b(m), where -N/2 <_ n, m _<_ N/2 - 1. 
Truncate these Fourier representations, keeping only the central two- 
thirds of the spectrum. 

J a(n) , |n| < N/3 

\ 0 , | n ] > N/3 


a(n) 


% 


(m) 


b(m) , |m| _< N/3 

0 , |m | > N/3 


(3-9) 


(2) Compute a(x ) and b(x ) from the truncated Fourier repre- 

J J 

sentation. 

(3) Form the product a(x ) b(x ) in physical space. 

^ j 

(4) Take the Fourier transform of the result, keeping only the 
central two-thirds of the spectrum. 


With this method of removing alias error, one-third of the mesh 
points in each direction are wasted in the one-dimensional case. In the 
three-dimensional case, 19/27 of the total mesh points are wasted. Con- 
sequently, this method is costly to use, and an alternative method of 
removing alias error will be described next. 
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C. Alias Removal By Phase Shifting 
We can rewrite equation (3-8) as 


C(s) 


Tj a(n) b(m) + ^ a(n) b(m) 

rrhn=s n+m*s±N 


(3-10) 


The first term on the right-hand side is the alias-free part of the 
result, while the second term on the right-hand side contains the terms 
that are aliased to wave-number s. 

The shift theorem tells us that if F(x) has the Fourier trans- 

a <Jrh «• 

form F(k), then F(x-h) has the Fourier transform e F(k) . This 
can be used as the basis of a method of removing aliasing error. The 
basic steps of phase shif ting method are : 


(1) Given a(n) and b(m) , form a(x ) and b(x ) on the 
shifted mesh. 


N/2-1 

E 

n— N/2 


Wxj) - £ Mn> . lnh(J+l> .SUj) - £ Mm) « 


N/2-1 

E 

mp-N/2 


imh(j+l) 


(3-11) 


where h 
( 2 ) 
(3) 

transform 


is the mesh size. 



Shift the result back to the original mesh. The Fourier 
of resulting function is: 


C’(s) 


i.h| 


E 


a(n)e 


inh 


b(m)e 


imh 


|n+m=s 


Y] a(n) e 
nfm“8±N 


inh 


b(m)e 


imh i 


“ 2 a(n)b(m) + e ±iNh a(n)b(m) 

n+m-s n+m»s±N (3-12) 


(4) Repeat steps (1), (2), and (3) with h replaced by h' * h + 
u/N. Equation (3-12) becomes 

C”(s) « a(n)b(m) + a(n)b(m) 

n4-m=s nrHn=*s±N 


* ^ a(n)b(m) - e iiNh ^ a(n)b(m) 

n+m=s n+m=s±N 


(3-13) 
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(5) The alias-free result is obtained by combining (3-12) and 
(3-13) 

^ a(n)b(m) « ~ (C'(s) + e"(s)J (3-14) 

n+m=s 

Ihe advantage of this method is that it does not waste any mesh 
points as in the ”2/3 rule". It doubles the computation but is far more 
efficient . 

3.3 Time Advancement 

In full-turbulence simulations we are looking for time-accurate 
solutions to the equations of motion. An explicit method of advancing 
the solution in time is a natural choice. In particular, we used a 
second-order Runge-Kutta method to advance the solution in time. The 
algorithm to advance the solution of du/dt ” f(u) from the time step 
n to step n + 1 is 

u * - u ( n > + A tf{u* n) J (3-15) 


(n+1) 

u 



(3-16) 


where At is the time step and * denotes the intermediate step. 

The accuracy of this method can be found by introducing the 
representative equation du/dt ■ Xu. Replacing f(u) by Xu in (3-15) 
and (3-16) and substituting (3-15) into (3-16), we obtain 

JnfD _ Q + XAt + j (XAt)^jj u (n) (3-17) 

Equation (3-17) shows this algorithm is indeed second order accurate. 

The stability limit of this method is the locus of all x 
satisfying the following relationship 

|| (XAt) 2 + (XAt) + l| - |e l9 | (3-18) 
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where 


Rewriting (3-18) 


0 < 0 < 2ir 


j (XAt) 2 + (XAt) + [1 - e i9 J - 0 , 0 < 0 < 2v (3-19) 

The stability diagram for this method is shown on Figure 3-1 . The time 
step is variable and is determined by setting the Courant number to be 
0.05. 


3.4 Initial Conditions 

There are several ways to generate the initial velocity field (Bar- 
dina et al., 1983; Shiranl et al., 1981; and Rogallo, 1981). We follow 
Rogallo’s method to produce an Initial velocity field which is nearly 
isotropic in wave space, divergence-free, and has the desired spectrum. 

None of the initialization processes mentioned above produces the 
higher-order velocity statistics of real turbulence. We therefore let 
the initial field decay until the higher order statistics and the proper 
energy cascade were developed before introducing compression. 

The Isotropic fields produced after the initial development period 
were used as the starting conditions for various straining runs. 
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Fig. 3-1. The stability diagram for second- 
order Runge-Kutta method. 
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Chapter IV 

COMPUTER PROGRAM DEVELOPMENT ON THE CYBER 205 


In this chapter, the CYBER 205 computer and the program used to do 
the simulations are briefly described. The data-management strategy 
chosen, which takes full advantage of this computer, is presented. The 
performance of the computer code and the validation checks made are 
described . 

4.1 The CYBER 205 Computer 

The simulations were done on two different Control Data Corporation 
CYBER 205 computers, one at Colorado State University and a larger one 
with twice the memory and number of processor units at Control Data Cor- 
poration, Minneapolis. The CSU computer has two million 64-bit words of 
central memory. Its central processor unit (CPU) contains a scalar pro- 
cessor and a vector processor with two pipelines. 

Figure 4-1 shows the performance of a two-pipe machine for addition 
and multiplication as a function of vector length. The asymptotic 
performance, which requires a vector length of 65535, is 100 million 
floating-point operations per second (Mflops) for 64-bit arithmetic and 
200 Mflops for 32-bit arithmetic. 

It is obvious that the performance improves with vector length. 
Vector length 1000 (64-bit case) or 2000 (32-bit case) is required to 
reach 90% of the asymptotic performance. Constructing a code which uses 
long vectors is therefore important if maximum performance is to be 
obtained from the machine. 

4.2 Data Management 

Based on the "longer vector gives better performance” philosophy, 
we chose to do the discrete Fourier transforms in parallel. This will 
be explained in detail later (also see Wu et al. 1983). 

In Figure 4-2, NX, NY, and NZ are the number of mesh points of 
the computational box in the x, y, and z directions, respectively; 
MY and MZ are called "pencil size". 
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On the first sweep, HZ x-y planes of data are Fourier transformed 
in the y direction; these are done in parallel. Ihe transform length is 
NY, but by doing them in parallel, a vector length of 3 x NX/2 x HZ is 
achieved; the factor 3 is due to the simultaneous processing of three 
velocity components, while the factor 1/2 arises because only half of 
the wave space modes are needed in order to represent a real function. 
To accomplish the vector izat ion, it is useful to lump every dependent 
variable into a single big array. The main array is called DATA (NX/2, 
NY, NZ, 4, 2); the dimensions represent x, y, z, a dependent variable 
index, and real and imaginary parts of a complex number. 

On the second sweep, MY x-z planes are processed. Fourier trans- 
forms in the z and x directions are done on this sweep. The vector 
lengths are NX/2 x MY x 3 and NZ x MY x 3, respectively. 

A CYBER 205 vector is defined as a contiguous set of memory loca- 
tions. Since the two sweeps are in different directions, an array 
transpose has to be done between sweeps and within the second sweep in 
order to keep processed data in a contiguous set of memory locations. 
The transpose is done by using gather instructions. The gather instruc- 
tion puts array elements which are at various locations into a contig- 
uous set of memory locations. An index vector is needed to pick up the 
desired elements. The Q8VGATHR function (64-bit) or the Q8VXT0V 
subroutine (32-bit), provided by the manufacturer, is used to do the 
transposing. As the array gets bigger, so does the index vector 
length, and an appreciable amount of overhead working space is needed. 
In the 64 x 64 x 64 run with pencil size 32 x 16, the index vector 
has 17,408 elements. 

4.3 Computer Performance 

The performance data based on a count of the number of operations 
per time step are presented in Table 4.1. The mesh size is given in 
column 2 (each node requires seven words of data storage). Hie pencil 
size is given in column 3; this, together with mesh size, determines the 
vector length shown in column 4. The computational precision is given 
in column 5, the CPU time in column 6, and the CPU computation rate in 
column 7. The I/O time per time step in seconds is meaningful only for 
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Table 4.1 
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runs with virtual memory paging. Explicit I/O would reduce I/O time 
considerably, but we did not use it. 

Comparing Runs 3 and 4, and Runs 5 and 6 in Table 4.1, it is found 
that the CPU time for a 32-bit (half) precision run is 60% of that for 
the corresponding 64-bit (full) precision run. We kept track of the 
timing in the transpose part of the code and found an interesting fact. 
In full-precision runs, the transpose takes 15% of the CPU time; 85% of 
the CPU time is spent in floating point operations. In half-precision 
runs, due to the lack of a half-precision gather utility, the transpose 
takes the same time as in full-precision runs, while the floating-point 
operations require only half of the full-precision CPU time. Conse- 
quently, for the half-precision run, the transpose takes 25% of the 
total time. 

Detailed timing from Run 8 shows that 51% of the CPU time is spent 
in the FFT subroutine, which contains 78% of the floating point opera- 
tions. In other words, the FFT operates at 157.6 Mflops. The remaining 
22% of the floating-point operations are executed at 95 Mflops, due 
mainly to shorter vector lengths and the occurrence of IF statements. 

Table 4.2 shows the performance data obtained from a CYBER 205 with 
four vector pipelines. The computation rate of a four-pipe machine is 
twice as fast as that of a two-pipe machine. Mote that Run 13 is only 
1.7 times faster than Run 8. This is again due to the transpose taking 
the same amount of time on both machines. 

In summary, the present code is fully vectorized to take maximum 
advantage of the capabilities of the CYBER 205 computer. It can handle 
up to 64 x 64 x 64 mesh points in core while running at 100 Mflops on 
a two-pipe machine and 190 Mflops on a four-pipe machine. This compu- 
tation rate is among the fastest that can be obtained from present 
supercomputers. To obtain it, the vector length must be long enough and 
half-precision computation must be used. Hie slowness of taking the 
transpose is the key factor in preventing one from getting better per- 
formance from the computer. The slow I/O transfer rate prevented us 
from using a 128 x 128 x 128 mesh. 
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4.4 Tests of the Code 


Several tests were performed in order to qualify the code. Two 
important checks are described below: 


A. Decay of Isotropic Turbulence 

A 64 x 64 x 64 run was performed with the mean strain rate tensor 

U set equal to zero. Ihe three-dimensional energy spectrum was ob- 
* »J 

tained by summing the energy in spherical shells of inner radius k and 
outer radius k + Ak, where k is the magnitude of the wavevector and 
Ak is the difference between the magnitudes of the two nearest neighbor 
wavenumbers. Ihe initial energy spectrum was a top-hat spectrum. The 
flow evolved from the artificial initial spectrum to a realistic low 
Reynolds number spectrum as the simulation proceeded. As can be seen 
from Figure 4-3, the normal Reynolds stresses are slightly anisotropic 
at low wavenumbers and isotropic at high wavenumbers. This is due to 
the small number of modes at low wavenumbers. The 3-D dissipation 
spectrum spectrum D(k) is defined as 


D(k) - J±k 2 E(k) (4-1) 

P 

and its components D (k) « k E(k) are shown in Figure 4-4 • The 

eta P 

peak of the dissipation spectrum is located well inside the resolvable 
wavenumber range. This means that the dissipation scale of the flow is 
resolved . 

The longitudinal and lateral one-dimensional spectra of the veloc- 
ity field are shown in Figure 4-5. The one-dimensional spectrum of the 
velocity field in the k^ direction is defined as 



(kj ) 



*(k) dk 2 dk^ 


(4-2) 


As shown, the spectra drop several orders of magnitude frow low wave- 
number to high wavenumber. This again demonstrates that the dissipation 
scales are well resolved. That E 22 (kj) - shows that the vel- 
ocity field is isotropic. 
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The two-point correlations and the one-dimensional spectra are 
Fourier transform pairs. They are defined as 


< u ! (x) u' (x + j) > 

Q (y) • ' ' — 

ij - < u^(x) u*(x) > 


(4-3) 


The close agreement between the two lateral correlations, 022^ x l*®*^) 
and Qj^Cxj , 0,0) , in Figure 4-6, shows that the flow is isotropic. 


The time evolution of the turbulent kinetic energy, dissipation 
rate, integral length scales, and Taylor micro-scales are all in good 
agreement with both experiments (Comte-Bellot and Corrsin, 1971) and 
previous simulations (Shiran! et al., 1981). These results gave confi- 
dence that the code is capable of simulating homogeneous isotropic tur- 
bulence. The particular flow field shown in this section was also used 
as the Initial condition for the three-dimensional straining runs. 


B. Plane Strain 

The second test of the code is isotropic turbulence subjected to 
uniform two-dimensional strain; one pair of sides of the computational 
box moves inward at the same rate a second pair moves outward, while the 
third pair remains stationary. This type of flow has been investigated 
experimentally by Townsend (1954) and by Tucker and Reynolds (1968); and 
numerically by Rogallo (1981). 

Three different strain rates and two different choices of strain 
axes were tested, namely: compression in the X 2 direction and stretch 
in the Xg direction, and compression in the x^ direction and stretch 
in the X 2 direction. The time development of the components of the 
turbulent kinetic energy and structure parameters of a typical run are 
shown in Figures 4-7 and 4-8. Figure 4-9 shows the time development of 
the components of the turbulent kinetic energy for two different choices 
of strain axes. 

The simulated results are in good agreement with Rogallo' s results. 
The differences can be attributed to different turbulence Reynolds num- 
bers. The close agreement between two different choices of strain axes 
shown in Fig. 4-9 demonstrates that the computer code does not have any 
direction preference. 
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From the results of the two tests described In this section, we 
have confidence that the computer code Is operating properly and ready 
to use to simulate compressed turbulence. 

4.5 General Description of the Simulations 

As described in Section 3.4, the simulation starts with an arti- 
ficial initial field. Therefore, there is an initial period in which 
the simulated flow fields cannot be treated as true turbulence. During 
this period, an energy cascade is gradually developed; the spectra be- 
come realistic, and higher velocity statistics reach asymptotic values. 

After the initial "relaxation" period, the velocity field is stored 
for later use as the initial condition for straining runs. During the 
straining process, the simulated flow fields can be regarded as true 
turbulent flows. Statistics are extracted from these flow fields. 
These will be discussed extensively in Chapters V and VI. This simula- 
tion is stopped when the scales of motion grow too large for the compu- 
tational box. 

Beyond this point, the flow fields do not accurately represent true 
turbulence, because the eddies are influenced by the imposed periodic 
boundary conditions. The region of validity of the simulation can also 
be monitored by examining the dissipation spectra. If the peak of the 
spectra moves beyond the resolvable wavenumber range or energy accumu- 
lates at high wavenumbers, the small scales are no longer being accu- 
rately resolved. When resolution of either the large or small scales is 
lost, the simulation is stopped. 

In the next chapter we discuss the simulated results for isotropic- 
ally compressed turbulence. 
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VECTOR LENGTH 

Fig. 4-1 . The performance of a two-pipe CYBER 
205 for addition and multiplication. 



Fig. 4-2. Sketch of the computational box. 


39 






40 


Pig. 4-3. Three-dimensional energy spectra for pl Three-dimensional dissipation spectra 

isotropic-decay turbulence. * f or isotropic-decay turbulence. 
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Fig. 4-9. The evolution of the component energy 
ratio for two different choices of 
strain axes. 
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Chapter V 


HOMOGENEOUS ISOTROPIC TURBULENCE UNDERGOING ISOTROPIC COMPRESSION 

In this chapter we present results for homogeneous isotropic turbu- 
lence undergoing isotropic compression. Numerical simulation results 
are compared with the predictions of rapid distortion theory. Results 
for four different compression rates are presented and discussed. 

5.1 Rapid Distortion Theory 

Rapid distortion theory (RDT) is a linear theory which describes 
the response of turbulence to a rapidly applied mean field. This theory 
was originated by Taylor (1935) and was extensively developed by Batch- 
elor and Proudman (1954). It is valid only in the limit of extremely 
rapid distortions. A measure of rapidity is the ratio of the turbulence 
time scale (q /e) to the mean flow time scale (1/S). RDT is valid 
when 

Sq 2 / e » 1 


where S, q, and e are the mean strain rate, turbulent velocity, and 
dissipation rate respectively. 

When turbulence is distorted by a rapidly applied mean flow, i.e., 
the time scale of the applied strain is much shorter than the turbulence 
time scale, nonlinear turbulence interactions are unimportant, so the 
turbulence energy cascade cannot come to equilibrium with the applied 
strain. Also, the dissipation process is too slow to be Important dur- 
ing the straining. Under these conditions the nonlinear Interaction and 
viscous dissipation terms in the momentum equations may be neglected. 
Thus, RDT is a linear inviscid analysis of turbulence. 

RDT analysis discussed in this section follows Reynolds' work 
(1984); we are grateful for permission to use his results. In RDT, it 
is more convenient to work with vorticity field than the velocity 
field. The vorticity is defined as 

“i 2 e ijk u k, j (5_1) 


preceding page BLANK NOT RIMED 
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where uj, and u k are the vorticity vector, the completely 
antisymmetric tensor of rank three (the Levi-Civita symbol), and the 
velocity vector, respectively. We can relate the velocity field to the 
vorticity field by taking the curl of Eq. (5-1). 


e ,o>, m u - u 

pqi i,q q»pq P,qq 


(5-2) 


The dynamic equations for the vorticity can be derived by taking 
the curl of the momentum equations : 


d(i>£ 

— + u iHt5 - S ±i - Wi S kk + £ (5-3) 

where S„ * | [u^ j + u j and p and y are functions of time 

only. The first two terms on the right-hand side of Eq. (5-3) are the 
vortex stretching terms and the last term accounts for viscous diffu- 
sion. 

All flow quantities in Eq. (5-3) can be decomposed into mean and 
fluctuating parts, that is 


“i “ W i + W i 


(5-4) 


u. 


"i + „; 


Substituting Eq.(5-4) into (5-3) and taking the ensemble average, we get 

aw. 


i - 
+ U W 

at 3 i»j 


Y« - V* + 


(5-5) 


The equations for the fluctuating vorticity (Eq. (5-6)) are obtained by 
subtracting Eq. (5-5) from (5-3). Note that, in deriving Eq. (5-6), we 
have assumed (1) that nonlinear interactions and viscous dissipation may 
be neglected, (2) that the mean strain is irrotational and (3) that the 
turbulence is homogeneous, these assumptions are applicable in the iso- 
tropic compression case. 


aw’ 

at" 


+ u w 


3 i,j 


w’S 


3 ij 


— w* S 
i kk 


(5-6) 


46 



Eq. (5-6) can be transformed to Lagrangian coordinates using the trans- 
formation of Chapter II. In the Isotropic compression case, Eqs. (2-21) 
and (2-34) are used to reduce Eq. (5-6) to: 



at 


2V 

-P-w' 

x (t) i 
P 


(5-7) 


where Vp is the compression speed and Xp(t) is the box length at 
time t. Ihe solution of Eq. (5-7) is 

/ x ( 

wj(x' ,t) « w^(x’ .0) / 

where L Q is the Initial box length, and Xp(t)/L Q is the instantane 1 
ous total strain ratio. 


r 


(5-8) 


The turbulent velocity field is found from the turbulent vorticity 
field via Eq.(5-2). 


u i,kk " " e ikj w j,k 

where continuity constraint u i i “ ® * s applied. 
Eq.(5-9) to Lagrangian coordinates, we get 


(5-9) 

Transforming 



(5-10) 


Equations (5-10) are solved by using Fourier transforms. Taking the 
Fourier transform of Eq.(5-10), we get. 




+ k’ 2 + k^ 2 j 


( 


k£w’(k',tj - k^w’^k' ,tj 
k’w’(k*,tj - k£w*(k’,tj 



(5-11) 
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where k' is the Fourier wave vector corresponding to x’ . We define 
Ej^Jk* ,tjand ^Jk' ,tj as the spectra of u’u* and w'w' , respec- 

J J A J * J 

tively. 






(5-12) 


** A 




i 


(5-13) 


where * denotes the complex conjugate. From Eqs. (5-11), (5-12), and 
(5-13), E^j^k* ,tj can be expressed as 

E llLlE'» t J " 2 ^,2 .,2 2 l k 3 *22 + k 2 *33 " 2k 2 k 3*23-* (5_14) 

l K l + k 2 + k 3 J 

and from Eqs .(5-8) and (5-13) ,tj can be expressed as 

x " 4 

♦ ijli'.tj ' ♦ylt’-OJ < 5 - 15 ) 


*22^’ »0) is obtained from Eqs. (5-13), and (5-12) 

♦22 * k 3 2 EnUc',0j + kJ 2 E 33 Qt’,Oj - 2kJk’E 13 yc’,0j (5-16) 


If the initial turbulence is isotropic, then 




(5-17) 


Combining Eqs. (5-14), (5-15), (5-16), and (5-17) we obtain 

-2 f 2 . .2 . .2 . .2 . .2 . 2 


V - ,tJ 


w ,0) ( l p (t) \ 

4»k' 2 \ L o / 


2 2 


_k; |k* -k; j + k; ^ -k; j + 2k* k* 


] 


, ,2 , ,2 ,,2 2 
i k ; + k 2 + k 3 J 

(5-18) 

The spectra of E 22 and E33 can be found by permuting the indices. 


Equation (5-18) is the most important result of RDT analysis for 
isotropically compressed turbulence. It provides the spectrum of the 
turbulence in terms of the initial spectrum and the Instantaneous total 
strain ratio. Any other quantity dependent upon the spectrum can be 
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deduced from Eq.(5-18). For example, the evolution of the three- 
dimensional dissipation spectrum D(k' ,t) is derived as follows: 

D(k' ,t) is defined as in Equation (4-1), 


D(k, ,t) - k ,2 (t) E(k',t) 


Substituting Eqs. (2-21), (5-19) and k'(t) * 
(4-1), we get 

x \-0.9-(-3)+(-2)+(-2) 


\“1 




k'(0) into Eq. 


D(k* ,t) 


{■<) 

/ x V 

l L / 

\ o' 


k’ 2 (0) E(k’,0) 


x V-1,9 


D(k’ ,0) 


Other important results from RDT analysis for isotropic compressed tur- 
bulence are summarized below: 


Three-dimensional energy spectrum: E(k',t) 


-2 

(4) 


E(k' ,0) 
-1,9 


Three-dimensional dissipation spectrum: D(k’ ,t) * 


One-dimensional energy spectrum: Ej^k^,tj 


(?) 

-2 

(4) V k i’°J 


D(k' ,0) 


Two-point auto-correlation: Q^x^.tj * Q^x^.Oj 

-2 

Turbulent kinetic energy: k(t) * (l^) ^(0) 


(5-19) 


Dissipation rate: e(t) 


(4) 


-1.9 


:( 0 ) 


v 2.05 


Kolmogorov length scale: nCb) “ ( XT') n(0) 


Taylor microscales: x^(t) 


t* *ii<°> 
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Integral length scales; L,. (t) 

ij ,m 

Turbulent Reynolds number; Re (t) 

£ 


x 

L 

o 


L ij >m 


( 0 ) 





5.2 Simulation Results 

Four simulations of isotropically compressed turbulence were per- 
formed. All started with the same isotropic turbulence initial condi- 
tion. but the compression rate differed in each case. Some details of 
the isotropic turbulence initial condition are shown in Table 5.1. 
Three independent dimensionless parameters can be formed from the time 
t , the strain rate S( t), and the turbulence velocity and time scales, q 
and q^/e» respectively. These are taken to be 

exp S(t’) dt'^ , Sq /e, q /ev* 

which are the total strain ratio, the ratio of the turbulence and strain 
time scales, and a turbulence Reynolds number, respectively. Since 
homogeneous turbulence is a time-developing flow, these dimensionless 
parameters change during a simulation. Table 5.2 shows the range of 
these parameters covered in each simulation. 

Run SQF had the fastest compression rate among these four cases. 
The ratio of turbulence and strain time scales is sufficiently high that 
the rapid distortion limit is applicable. The results of run SQF agree 
with RDT analysis (Eq .<5—19) ) ; this is demonstrated below. 


A. Description of the Rapid Distortion Run 

In this section, we present results obtained from the rapid distor- 
tion run (SQF). The mean strain-rate tensor is 

( S(t) 0 0 

0 S(t) 0 

0 0 S(t) 

where S(t) “ V /X ( t ) . The results are shown in Figures 5-1 through 

r r 

5-11 and discussed below. 
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Table 5.1 


Some Details of the Initial Condition 
of the Isotropically Compressed Turbulence Simulations 
(consistent units) 


Turbulence Velocity q 

0.2856 

Dissipation Rate e 

0.0324 

Integral Length Scale Ljj ^ 

0.9791 

Taylor Microscale Aji 

0.3444 

Kolmogorov Length Scale q 

0.07454 

Kinematic Viscosity v 

0.01 
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Figures 5-1 through 5-4 show the evolution of the lateral and long- 
itudinal one-dimensional spectra of the velocity field; these are de- 
fined by Eq.(4-2). The initial lateral and longitudinal spectra are 
shown in Figures 5-1 and 5-3, respectively. The spectra obtained from 
the simulation and RDT analysis when the total strain ratio is 0.552, 
are shown in Figures 5-2 and 5-4, respectively. The agreement is excel- 
lent. 

Figures 5-5 and 5-6 show the longitudinal and lateral two-point 
auto-correlations at the initial time, and when the total strain ratio 
is 0.779 and 0.552. They are given in the Lagrangian coordinates. 
These figures show that, In accord with RDT analysis (Eq.(5-19)), the 
two-point auto-correlation in Lagrangian coordinates does not change 
during the simulation. 

Figures 5-7 and 5-8 show the time evolution of twice the turbulent 
kinetic energy (q^) and the dissipation rate (e) from the simulation 
and RDT analysis, Again, the agreement is excellent. 

The time evolution of the integral length scale, Taylor microscale, 
and Kolmogorov length scale from the simulation and RDT analysis are 
shown in Figs. 5-9, 5-10, and 5-11. The integral length scale, which is 
the integral of the two-point velocity auto-correlation, shrinks line- 
arly with time in Cartesian coordinates. So does the Taylor microscale. 
However, the Kolmogorov length scale decreases much faster than the in- 
tegral length scale. This is due to the rapid decay of the kinematic 
viscosity, which is, in turn, due to the increasing density and temper- 
ature during compression* This results in the small length scale 
decreasing very rapidly during the compression and is important in 
understanding some of the results presented later. 

In summary, the agreement between the results of RDT analysis and 
computer run SQF is excellent. This demonstrates that run SQF indeed is 
a simulation of rapid distortion and, taken together with the results of 
the isotropic decay run, that the code is performing correctly. 
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B. Summary of Isotropic Compression Simulations 

Important data from all four isotropic compression simulations are 

compiled in Table 5.2 and discussed in this section. These simulations 

cover a wide range of initial values of |s|q /e, from 0.1 to 47.04. 

As shown above, the fastest compression simulation agrees with rapid 

distortion analysis. On the other hand, for extremely slow compression 
2 

(Js|q /e ♦ 0j, the effect of strain becomes negligible and the flow 
behaves like the isotropic turbulence decay of Chapter IV. 

In the region between these two extremes, there are no known ana- 
lytical solutions. The effect of |s|q /e on the various turbulence 
statistics in this region is discussed below. During the compression 
stroke of an internal combustion engine, the value of |s|q /e is in 
the range 0.05-0.5. So, runs SQH and SQl are in a range applicable to 
engines . 

Figures 5-12 and 5-13 show the behavior of the Reynolds number and 
time scale ratio as functions of the total strain, which is an approp- 
riate nondimens ional time for this flow, for the four simulations. The 
four simulations have the same initial condition shown in Table 5.1 but 
different compression speeds. Both the Reynolds number and time scale 
ratio increase with time. For cases with small initial |s|q /e* the 
rate of increase of the Reynolds number is less than for cases with 
larger initial time scale ratio. As can be seen from Fig. 5-13, a wide 
range of |s|q /e has been covered by these four simulations. 

Figure 5-14 shows the three-dimensional energy spectrum at total 
strain ratio 0.785 in run SQH. The 3-D energy spectrum evolves as the 
simulation proceeds. No inertial subrange (a region of the spectrum 
with slope -5/3) appears in the 3-D energy spectrum because the Reynolds 
number is low. The flow field remains isotropic throughout the simula- 
tion. The corresponding spectrum from the rapid distortion simulation 
is plotted for comparison. The difference between these two spectra can 
be attributed to the effects of viscous dissipation and nonlinear inter- 
actions. 

Three-dimensional dissipation spectra for runs SQH and SQF at total 
strain ratio 0.785 are shown in Fig. 5-15. The peaks of the spectra are 
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located well inside the resolvable wavenumber range. This means that 
the dissipation scale of the flow is resolved. The difference between 
these two spectra can be attributed to the fact that run SQH is far from 
rapid distortion limit. 

The lateral and longitudinal two-point correlations of runs SQH and 
SQF at total strain ratio 0.785 are shown in Lagrangian coordinates, in 
Figs. 5-16 and 5-17. As can be seen, curves from run SQH are fatter 
than those from the rapid distortion simulation (run SQF). Conse- 
quently, the integral scales from this simulation are bigger than those 
from the rapid distortion simulation. In other words, nonlinear and 
viscous effects cause the integral scales to decrease less rapidly than 
in the rapid distortion case. 

Figures 5-18 and 5-19 show the history of the turbulent kinetic 
energy and its dissipation rate. When the compression rate is large, as 
in runs SQF and SQG, the flow is immediately affected by the mean strain 
and both quantities increase throughout the simulation. If the initial 
value of |s|q /e is 0.1 or less, the mean strain is not strong enough 
to immediately alter the structure of the flow and the flow retains the 
character of decaying isotropic turbulence during the initial period. 
After the strain has been applied for a longer period, the flow grad- 
ually reflects the effect of mean strain. There is a transition region 
in which the decay and strain effects are approximately balanced. In 
all cases, when the strain has been applied for a long enough period, 
both the turbulent kinetic energy and its dissipation rate increase. 

Figures 5-20 through 5-22 show the evolution of integral length 
scale, Taylor microscale, and Kolmogorov length scale, respectively. 
These length scales all behave in roughly the same way. If the com- 
pression is fast enough, the length scale immediately begins to 
decrease; if the mean strain is weaker, the length scale may grow 
initially as in the case of decaying turbulence. In the latter cases, 
after the strain has been applied for a while, the length scales start 
decreasing. Note that the Kolmogorov length scale decreases most rap- 
idly. The increasing density is responsible for the dramatic decrease 
of the kinematic viscosity. In most cases, the simulations had to be 
stopped because the small scales could no longer be resolved. 
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Fig. 5-4. The longitudinal one-dimensional 
Fig. 5-3. The initial longitudinal one-dimen- energy spectra of run SQF and RDT 

sional energy spectrum of run SQF. analysis at total strain 0.552. 
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Fig. 5-7. The evolution of twice the turbulent Fig. 5-8. The evolution of the dissipation rate 

kinetic energy of run SQF and ROT of run SQF and RDT analysis, 

analysis. 
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Pig. 5-9. The evolution of a longitudinal inte- Fig. 5-10. The evolution of a Taylor microscale 
gral length scale of run SQF and RDT of run SQF and RDT analysis, 

analysis. 
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Fig. 5-11 . The evolution of the Kolmogorov Fig. 5-12. The evolution of the turbulent 

length scale of run SQF and RDT Reynolds number for four isotropic 

analysis. compression simulations. 
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Fig. 5—19 . The evolution of the dissipation | Fig. 5-20 . The evolution of a longitudinal In- 
rate for four Isotropic compression tegral length scale for four Iso- 

slmulatlons. tropic compression simulations. 
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Pig.. 5-21. The evolution of a Taylor microscale Fig. 5-22. The evolution of the Kolmogorov 
for four Isotropic compression length scale for four isotropic 

simulations. compression simulations. 



Chapter VI 


RESULTS FOR HOMOGENEOUS TURBULENCE UNDERGOING 
QNE-DIMENSIQNAL COMPRESSION 

In this chapter we present the results for homogeneous turbulence 
undergoing one-dimensional compression and expansional, irrotational , 
axisymmetric strain. Four one-dimensional compression and three axisym- 
metric strain simulations are presented and discussed in some detail. 
The axisymmetric strain simulations were carried out by Lee and Reynolds 
(1985); we are grateful for permission to use their results prior to 
publication. 


Description of The Flows 


In this section we describe the flows that are simulated numeric- 
ally. We consider a fixed mass of turbulent fluid contained within a 
rectangular parallelepiped. If one opposed pair of sides of the paral- 
lelepiped moves inward, as illustrated in Figure 2-2, the fluid in the 
box undergoes a one-dimensional compression similar to that experienced 
in an engine cylinder with a flat piston. If one opposed pair of sides 
moves inward at twice the rate the other two pairs move outward , the 
fluid experiences an expansion type of axisymmetric strain. This case 
is illustrated in Figure 6-1. Both flows are axisymmetric. Indeed, the 
one-dimensional compression can be considered as a combination of the 
axisymmetric strain and the isotropic compression studied in the preced- 
ing chapter. Formally, one can decompose the strain-rate tensor for the 
one-dimensional case into those for the other two cases : 



Tables 6.1 and 6.2 present the range of independent dimensionless 
parameters covered in each simulation. The run identification is given 
in Goluran 1. The initial and final |s|qVe are given in Columns 2 and 
3, respectively. A wide range of |s|q^/e was covered in each type of 
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simulation so that the effect of strain rate on turbulence could be 
studied. Columns 4 and 5 show the initial and final turbulent Reynolds 
numbers for each simulation. Each run covers much the same range of 
Reynolds number. The anisotropy tensor b^j is defined as 


R 


'U 


“ij °ij 
R kk 3 


(6-2) 


where Rjj is the Reynolds stress tensor and 
delta. Three invariants can be formed from b 


°ij 


is the Kronecker 




’ii 


(6-3) 


-211 


b ij b ji 


(6-4) 


3IH " bijb jk b ki (6-5) 

The first invariant (Eq. (6-3)) is zero as a consequence of the defini- 
tion of bjj . II and III are measures of the anisotropy of Rjj or the 
magnitude of bjj . The initial and final values of II and III are shown 
in Columns 6 through 9 of Tables 6l.l and 6.2. Since both types of strain 
are axisymmetric , b£2 * b^j » -bjj/2. The final total strain ratio is 
shown in Column 10. 

The same initial state was used for every simulation and was taken 
from a simulation of decaying, homogeneous, isotropic turbulence that 
had been run long enough to allow the energy cascade to become estab- 
lished. The one-dimensional, compressed turbulence simulations were 
performed on a grid of 64 x 64 x 64 mesh points. Initially, the com- 
putational cell has * 2 ax£ ■ 2ax^, where is the compression 
axis. The total strain ratio achieved is approximately 0.25; higher 
strains result in Inadequate resolution, because the box becomes too 
small in the compression direction. The expansional, axisymmetric 
strain simulations were done by Lee and Reynolds (1985) and were per- 
formed on a grid containing 128 x 128 x 128 mesh points. At the ini- 
tial state, the computational cell has &x^ * 2 v2ax^ “ 2/2ax ^ . The 
total strain ratio achieved is close to 2. Results of these two types 
of simulations are described and discussed in the following sections. 
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6.2 Simulation Results 


A. Summary of One-Dimensional Compression Simulations 

Four simulations of one-dimensional, compressed turbulence were 
performed. All started with the same isotropic turbulence initial 
condition, but differed in compresssion rate. Some details of the 
initial field are shown in Table 6.3. Some overall data characterizing 
the four one-dimensional compression simulations are compiled in Table 
6.1 and discussed in this section. 

The mean strain-rate tensor for one-dimensional compression is 



where S(t) ■ Vp/ x p (t) , Vp is the compression speed, and Xp(t) l g 
the instantaneous box length. Summarized results are shown in Figures 
6-2 through 6-15 and discussed below. 

Figures 6-2 and 6-3 show the behavior of the Reynolds number and 
time-scale ratio as functions of the total strain, which is an approp- 
riate nondimenslonal time for this flow, for the four simulations . As 
can be seen in Figure 6-3, the time-scale ratio increases with time in 
each case. A wide range of |s|q^/e was covered by these four simula- 
tions to allow study of the effect of time-scale ratio on turbulence. 
As shown in Figure 6-2 for the cases with small initial |s|q^/e (runs 
ODD and ODE), the Reynolds number decreases in the initial stages of the 
simulation, and then increases. In the large Initial j S | q 2 / e cases 
(runs ODB and ODC), the Reynolds number increases monotonically with 
time. 

The evolution of the three-dimensional energy and dissipation spec- 
tra of a typical run are shown in Figures 6-4 and 6-5. These particular 
spectra are taken from run ODD at the initial time and at total strain 
0.494. The 3-D energy and dissipation spectra evolve as the simulation 
proceeds . Apparently, viscosity is dissipating the small scales while 
the large scales are absorbing energy from the applied strain. 

Figure 6-6 shows the three-dimensional spectra of the three compo- 
nents of the turbulent kinetic energy in run ODD at total strain ratio 
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0.494. Anisotropy occurs at the low and moderate wavenumbers. The 3-D 
energy spectum of Ejj is greater than that of the other two compo- 
nents, due to the fact that x^ is the compression axis. E 22 and E 33 
are almost equal, as expected. 

Figure 6-7 shows two 3-D energy spectra of runs ODB and ODD at 
total strain 0.847. Run ODB has a higher compression rate than run ODD, 
Both spectra evolve from the same initial spectrum, but differences 
develop. The differences can be attributed to the viscous dissipation 
of the small scales and nonlinear interactions. 

For runs ODB and ODD at total strain ratio 0.847, the longitudinal 
two-point correlations in Lagrangian coordinates are shown in Figures 
6-8 and 6-9. Note that x^ is the compression axis and X 2 is a non- 
compression axis. As can be seen, curves from run ODD are slightly 
broader than those from run ODB. Consequently, the integral scales from 
run ODD are larger than those from run ODB in all three directions. 
Nonlinear and viscous effects cause the integral scales to decrease less 
rapidly in the small | S | q^/e than in the large |s|q^/e case. 

Figures 6-10 and 6-11 show the history of twice the turbulent kin- 
o 

etic energy (q ) and its dissipation rate (e)* Both quantities show 
the same behavior as in isotropically compressed turbulence. When the 
compression rate is large, as in runs ODB and ODC, the flow is immedi- 
ately affected by the mean strain and both quantities increase through- 
out the simulation. If the initial value of |s|q^/e is 0.5 or less, as 
in runs ODD and ODE, the mean strain is not strong enough to Immediately 
alter the structure of the flow and the flow retains the character of 
decaying isotropic turbulence during the initial period. After the 
strain has been applied for a longer period, the flow gradually reflects 
the effect of mean strain and both quantities increase. 

Figure 6-12 shows the evolution of 1 >,j, which is defined by Eq. 
(6-2). It appears that there may be an asymptotic value of this quant- 
ity in each simulation, although the asymptote is not reached in the 
time covered by the simulations. The asymptotic value of bjj de- 
creases as |s|q^/e increases. Furthermore, for fixed total strain, 
b^j increases as time scale ratio ( |s|q^/e) decreases. The reasons 
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for these results, which appear to conflict with intuition, are not 
known at present. 

Figures 6-13 through 6-15 show the evolution of integral length 
scale, Taylor microscale, and Kolmogorov length scale, respectively. 
These length scales all behave in roughly the same way and are similar 
to those in isotropically compressed turbulence. If the compression is 
fast enough, the length scale immediately begins to decrease; if the 
mean strain is weaker, the length scale may grow Initially as in the 
case of decaying turbulence. In the latter case, after the strain has 
been applied for a while, the length scales start decreasing. Note that 
the rate of decrease of the Kolmogorov length scale in this flow is 
smaller than in the Isotropic compression cases, due to the fact that 
the density increases linearly with the Inverse of the total strain in 
this flow rather than cublcally, as in the isotropic compression cases. 

B. Summary of Incompressible Axisymmetric Expansion Simulations 

One-dimensional compression can be considered as a combination of 
isotropic compression and axisymmetric expansion (see Eq. (6-1)). If S 
is negative, the left-hand side of Eq. (6-1) is the mean strain-rate 
tensor for one-dimensional compression flow. The first term on the 
right-hand side is the mean strain-rate tensor for the isotropic com- 
pression, and the second term is the mean strain-rate tensor for axi- 
symmetric expansion. Although the mean strains are related by Eq. 
(6-1), the turbulence they generate may not be related to each other. 
The isotropic and one-dimensional compression simulations were described 
in detail in Chapter V and the preceding section. The relationship 
among the mean strains for the three flows suggests that there might be 
a relationship among the turbulence generated in the two compression 
flows and the one generated by the application of the mean strain cor- 
responding to the last term in Eq. (6-1). Tb study this, we need the 
data from the results of the axisymmetric expansion simulations of Lee 
and Reynolds (1985); they will be introduced in this section. The re- 
sults Introduced in this section will be used later in turbulence model 
testing. 
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The mean strain-rate tensor for axisymmetric expansion flow is 



where S is a positive constant. The turbulence is compressed in the Xj 
direction at twice the rate that it is stretched in the other two direc- 
tions. The mean velocity is divergence-free, so the density remains 
constant throughout the simulation. Figure 6-1 illustrates this type of 
flow. In Section 6.1 we present a general description of the simula- 
tions. Table 6.2 shows the range of independent, dimensionless param- 
eters covered in each run. Some details of the isotropic-turbulence 
initial condition are shown in Table 6.4. Summarized results are shown 
in Figures 6-16 through 6-23 and discussed below. 

Figure 6-16 shows the behavior of the Reynolds number as a function 
of the total strain ratio for each run. In the low-strain-rate run (run 
EXO) the Reynolds number decreases in the initial period of the simula- 
tion and later Increases. In the late stages of the simulations, the 
highest strain-rate run (run EXQ) has the lowest Reynolds number and the 
lowest strain-rate run has the highest Reynolds number. This contrasts 
with the 1-D compression flow. 

The behavior of the time-scale ratio (|s|q 2 /e) as a function of 
the total strain ratio is shown in Figure 6-17. As can be seen, a wide 
range of |s|q 2 /e is covered in these simulations. Since the mean 
strain rate (S) remains constant throughout each simulation, the 
change of the time-scale ratio is caused by the change of the turbulence 
time scale (q /e) only. It is found that the turbulence time scale 
increases during run EXO, which has the lowest strain rate. The turbu- 
lence time scale decreases slowly in the highest strain-rate run (run 
EXQ), and it is approximately a constant throughout run EXP. 

Figures 6-18 and 6-19 show the time evolution of twice the turbu- 
lent kinetic energy (q 2 ) and its dissipation rate (e). The behavior 
of both quantities is similar to that in one-dimensional compression 
cases. If the mean strain is strong enough, the structure of the flow 
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is immediately altered and both quantities increase throughout the 
simulation. Otherwise, the flow retains the character of decaying 
isotropic turbulence until it reflects the effect of mean strain. 

The evolution of bjj is shown in Figure 6-20. Its behavior is 
also very similar to that found in the one-dimensional compression case. 
The asymptotic value of bjj in run EXQ is 1/6, which agrees with RDT 
analysis. At the same total strain ratio, the flow which has the 

A 

smaller value of |s|q z /e is more anisotropic than the one that has 
larger |s|q 2 / e . 

Figures 6-21 through 6-23 show the evolution of integral length 
scale, Taylor microscale, and Kolmogorov length scale, respectively. 
The behaviors of these three length scales are similar to each other and 
to those of the one-dimensional compression flow. The behavior of the 
Kolmogorov length scale is closer to that of the other length scales in 
this flow (in contrast to the case of the one-dimensional compression 
flow), because the kinematic viscosity remains constant in this flow. 

Important results from three types of simulated flows — isotropic 
compression, one-dimensional compression, and axisymmetric expansion 
flows — have been given and discussed in this and the previous chapters. 
In the next chapter we shall use these results to test the validity of 
turbulence models. 
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Table 6.3 


Some Details of the Initial Condition 
of the One-Dimensional ly Compressed Turbulence Simulations 

(consistent units) 

Turbulence velocity q 0.9447 

Dissipation rate e 1 .5369 

Integral length scale 0.3227 

Taylor microscale Xjj 0.2146 

Kolmogorov length scale n 0.0385 

Kinematic viscosity v 0.0150 

Table 6.4 

Some Details of the Initial Condition 
of the Axisymmetrlc Expansion Flow Simulations 
(consistent units) 

Turbulence velocity q 0.4688 

Dissipation rate e 0.1931 

Integral length scale 0.2271 

Taylor microscale Xjj 0.1556 

Kolmogorov length scale n 0.0253 

Kinematic viscosity v 0.0043 
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Fig. 6-1 . Schematic description of homogeneous. 
Incompressible, axisymmetric expan- 
sion flow. 
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Fig. 6-6. The three-dimensional spectra of the Fig. 6-7. The three-dimensional energy spectra 

three normal energy components of run of runs ODB and ODD at total strain 

ODD at total strain 0.494. 0.847. 
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Pig. 6-8. The longitudinal two-point velocity Fig. 6-9. The longitudinal two-point velocity 

correlations in the compression correlations in the noncoopression 

direction of runs ODB and ODD at direction of runs ODB and ODD at 

total strain 0.847. total strain 0.847. 
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1/ (TOTAL STRAIN) 1/ (TOTAL STRAIN) 

Fig. 6-13. The evolution of the longitudinal 
, . _ integral length scale in the com- 

g. 6-12. The evolution of bjj for four one- press ion direction for four one- 

dimensional compression simulations. dimensional compression simulations. 
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TOTAL STRAIN TOTAL STRAIN 

Fig. 6-16. The evolution of the turbulent Rey- Fig. 6-17, The range of Jsfq^/e covered in 

nolds number for three axlsymmetric three axisymmetric expansional 

expansional strain simulations. strain simulations. 
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Pig. 6-18. The evolution of twice the turbulent Fig. 6-19. The evolution of the dissipation 

kinetic energy for three axisymmet- rate for three axisymmetric expan- 

ric expansional strain simulations. sional strain simulations. 
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Pig. 6-22. The evolution of a Taylor microscale Pig* 6-23. The evolution of the Kolmogorov 

for three ax i symmetric expanslonal length scale for three axlsymmetrlc 

strain simulations. expanslonal strain simulations. 




Chapter VII 

k- e TURBULENCE MODEL TESTING AND MODIFICATIONS 


Results from isotropic compression, axisymmetrie expansional 
straining, and one-dimensional compression simulations are used to test 
the validity of k-e turbulence models in this chapter. The model is 
found to be defective. To correct the deficiency, a three-equation 
(k-e-x) turbulence model which models four types of flows (isotropic 
decay, isotropic compression, axisymmetrie expansional straining, and 
one-dimensional compression flows) with a wide range of strain rates 
is proposed. A method of determining the model constants is developed 
and discussed. 

7.1 Introduction 

Since turbulence phenomena are highly nonlinear, they are not 
amenable to classical analytical approaches, and modeling needs to be 
used. As technology advances, greater accuracy is needed and more 
complex prediction methods, including more sophisticated turbulence 
models, are required. 

Methods of simulating turbulent flows can be classified according 
to the following scheme (Kline et al., 1981): 

1 . Correlations 

2. Integral methods 

3. One-point closure methods 

4. Ttoo- point closure methods 

5. Large-eddy simulation 

6. Full simulation 

As one moves down the list, each method requires less modeling than 
those above it. Hie range of flows that may be simulated with a single 
model broadens as the level increases. Turbulence models at each level 
can be validated using results from higher levels. Two-equation models 
(and k-e models in particular) based on one-point closure are among 
the most popular turbulence models at present. They are able to simu- 
late a wide variety of flows with reasonable accuracy, including flows 
that are difficult to treat with other one-point closure models. 
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The isotropic and one-dimensional compression simulations described 
in the preceding chapters are full simulations, i.e., numerical solu- 
tions of the exact Navier-Stokes equations. These results will be used 
to investigate both the physics of turbulence and turbulence models. In 
the next section we shall use the simulation results to test the valid- 
ity of k-e two-equation turbulence models; they will be shown to be 
inadequate. In Section 7.3, a new one-point-closure , three-equation 
turbulence model is proposed to overcome the shortcomings of k-e mod- 
els. The rationale for constructing the three-equation model is dis- 
cussed, a method for evaluating the model constants is described, and 
the performance of this three-equation model is presented. 

7.2 Testing k-e Turbulence Models 

A. Background 

The k-e two-equation model was developed by Jones and Launder 
(1972) and has been used in many applications. In addition to the 
equations for the mean flow, this model uses two partial differential 
equations which describe the evolution of the turbulence kinetic energy 
(k “ q^/2) and its dissipation rate (e)» In this model, the length 
scale is taken to be 

The equation which describes the turbulence kinetic energy can be 
derived by: taking the scalar product of Eq. (2-10) with the fluctuating 
velocity (u^) and taking the ensemble average. For homogeneous turbu- 
lence, the result is: 


dk 

■at " 


P - e 


(7-1) 


where P * - R TT is the rate of production of turbulence energy and 

ij i,j 

e • v u 1 u* is the rate of viscous dissipation of turbulent kin- 
i*j i»J 

etic energy; all quantities are per unit mass, 
is the Reynolds stress tensor. 


ij 


u’u' ; 

i i 


~ P R 


ij 


One can derive an equation for the dissipation rate (e) from the 
Navier-Stokes equations. The resulting equation is so complicated (El 
Tahry, 1983) that one has to model all of the terms in it. For homogen- 
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eous turbulence, the most commonly used model equation for the dissipa- 
tion rate (e) has the following form: 


'lF + C 3 eS kk 


where the are model constants. The first term on the right-hand 
side of Eq. (7-2) is the only one active in isotropic decaying turbu- 
lence, and the constant, C£, can be fixed by fitting it to that flow. 
The second term accounts for production by incompressible strain and/or 
shear, and Cj is adjusted to produce the correct behavior for incom- 
pressible homogeneous straining and shearing flows (Reynolds, 1976). 
The third term accounts for production due to dilatation; the evalution 
of C 3 requires data for compression flow; the results of Chapter V 
provide the necessary data (Reynolds, 1980). 

Watkins (1977) suggested that the constant C 3 be unity as a natu- 
ral outcome of his derivation. Hoult and Wong (1980) argued that the. 
angular momentum of the turbulence should be conserved in rapid compres- 
sion. Using this notion, Reynolds (1980) suggested that the the product 
of the turbulence length and velocity scales should remain fixed during 
fast compression, i.e.: 


constant 


(7-3) 


Tb obtain the proper behavior of turbulent length scale in rapid iso- 
tropic compression, Reynolds proposed C 3 ■ -2/3. 

Two sets of model constants are tabulated in Table 7.1. One was 
proposed by Launder and Spalding (1974) and extended by Watkins (1977) 
(hereafter denoted by LSW), and the other was proposed by Reynolds 
(1980). Note that not only are the magnitudes of model constant C 3 
different in the two sets, but the signs disagree. Resolution of this 
discrepancy was one of the original goals of this work. 
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Table 7.1 


The Values of the Model 

Constants 

in the k-e 

Models 


C 1 

C 2 

c 3 

C u 


Launder , 

Spalding, (LSW) 1.44 

Watkins 

1 .92 

1.00 

0.09 


Reynolds 1 .0 

11/6 

-2/3 

0.09 



To close equations (7-1) and (7-2), additional equations must be 
provided for Reynolds stress tensor R^j . In the k-e model, R^j is 
assumed to obey the Boussinesq constitutive equation: 


R ij " I k6 ij “ 2v tl S ij “4 S kk 6 ij ) 


(7-4) 


where 2k “ R^, *• y ^ ^j, v t is the eddy viscosity, 

and is the Kronecker delta. The eddy viscosity vt g* ven by: 


v 


t 



(7-5) 


where is chosen to be 0.09 to fit the ratio of shear stress to tur- 
bulence kinetic energy in local-equilibrium free shear layers, in which 
production and dissipation are equal. 

Tests of k-e models for isotropic compression, axisymmetric ex- 
pansion, and one-dimensional compression flows will be described below. 
In these tests, Eqs. (7.1) and (7.2) were solved with given initial 
values of k and e and the strain-rate tensor (Sjj) was prescribed 
as a function of time. The model (Eqs. (7-4) and (7-5)) for the Rey- 
nolds stress is used to compute the production. 


B. Isotropic Compression Flow 

The strain-rate tensor S.. for isotropic compression is 


’±5 


° ) 
S(t)/ 


S(t) 

0 

0 
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0 

S(t) 

0 


(7-6) 






where S < 0. The model (Eqs . (7-4) and (7-5)) for the Reynolds stress 
is not needed, because the turbulence remains isotropic throughout the 
simulation. Predictions of the two k-e models described above for 
four different strain rates, which correspond to runs SQF through SQI of 
Table 5.2, are shown in Figures 7-1 through 7-4, respectively. In the 
high strain-rate test (run SQF), both models predict k correctly but 
e is badly predicted. In this case, the production dominates the dis- 
sipation, so the prediction of energy is good because the production is 
correctly computed (the Boussinesq relation, Eq. (7-4), is not needed in 
this flow) . The dissipation is poorly predicted but has no effect on 
the prediction of the kinetic energy. In the moderate strain-rate runs 
shown in Figures 7-2 and 7-3, Reynolds' model predicts too high a dissi- 
pation rate causing the predicted energy (k) to be low. The LSW model 
predicts too low a dissipation rate and overpredicts the energy. In the 
slow compression rate test (run SQI), Reynolds' model produces good 
results while the LSW model overpredicts the energy. In these cases 
(except run SQF) , the energy production and dissipation are nearly 
equal, so accurate prediction of dissipation rate is required for a good 
simulation. 

The turbulence length scale is defined in k-e model as 

% - k 3/2 / e (7-7) 

and is supposed to represent the integral length scale. Figure 5-20 
shows the behavior of the Integral length scale in this flow. Figures 
7-5 and 7-6 show the behavior of the model length scale (t) predicted 
by the LSW and Reynolds' models, respectively. The LSW model predicts 
growth of the length scale during compression no matter how fast the 
flow is compressed. Morel and Mansour (1982) found that, in an engine 
simulation, the length scale predicted by the LSW model reached values 
several times the cylinder height near top dead center. This is clearly 
incorrect. The behavior of the length scale predicted by Reynolds' 
model is much more plausible. However, in order to get proper length 
scale behavior, Reynolds' model predicts too much dissipation in the 
fast compression case. 
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where S is a positive constant. The Reynolds stress model, Eqs. (7-4) 
and (7-5), is needed to model the production. Results for three differ- 
ent strain rates, which correspond to runs EXQ, EXP, and EXO of Table 
6.2, are shown in Figures 7-7 through 7—9, respectively. In the high 
and intermediate strain-rate tests (runs EXQ and EXP), predictions of k 
and e by both models behave badly. This is mainly due to overpre- 
diction of the production which is, in turn, due to overprediction of 
the Reynolds stress in the high and intermediate strain-rate cases. The 
Reynolds stress model, Eqs. (7-4) and (7—5), causes the production of 
energy to be proportional to the square of the instantaneous strain rate 
and predicts a sudden jump in Reynolds stress when strain is turned on; 
neither of these predictions is correct. In the low strain-rate test 
(run EXO), the predictions of both models are good and there is almost 
no difference between them. 

Figures 7-10 and 7-11 show the behavior of the turbulence length 
scale, defined by Eq. (7-7), predicted by the k-e models. Again, both 
k-e models predict growth of the length scale during the straining pro- 
cess no matter how large the strain rate. Figure 6-21 shows the evolu- 
tion of the integral scale in this flow. As can be seen, both models 
perform poorly. 


D. One-Dimensional Compression Flow 


The strain-rate tensor S., for one-dimensional compression is 



(7-9) 
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where S is negative. The Reynolds stress model (Eqs. (7-4) and (7-5)) 
is. used to calculate the production. Results for four strain rates, 
which correspond to runs ODB through ODE of Table 6.1, are shown in Fig- 
ures 7-12 through 7-15, respectively. In the high strain-rate test (run 
ODB) , both k and e are badly missed due to incorrect modeling of the 
Reynolds stress. In the intermediate strain rate case (run ODC), Rey- 
nolds' model predicts too much dissipation but the energy is in good 
agreement with the simulation results. The LSW model does not behave as 
well as Reynolds' model. In the slow strain-rate tests (runs ODD and 
ODE), Reynolds' model does an excellent job of predicting the turbulence 
energy (k) , dissipation rate (e)» and Reynolds stress (Rjj). The 
LSW model behaves better than in the high strain-rate tests, but is not 
as good as Reynolds' model. 

The evolution of the integral length scale for this flow is shown 
in Figure 6-13. Figures 7-16 and 7-17 show the behavior of the turbu- 
lence length scale predicted by the two models . The LSW model predicts 
growing length scales during compression which is not correct, while 
Reynolds' model shows plausible turbulence length-scale behavior (except 
in case ODB) . 

E. Comments on Turbulence Length Scale Model 

The turbulence length scale, which is supposed to represent the 
integral length scale, is modeled as k^^/e. Figures 7-18 through 7-20 
show the evolution of k^^/ e , based on the exact quantities obtained 
from the full simulations, as a function of total strain in isotropic 
compression, axisymmetric expansion, and one-dimensional compression 
flows, respectively. The longitudinal integral length scale in the com- 
pression direction for each flow is shown for comparison. As can be 
seen, this turbulence length-scale model behaves more or less correctly 
only in the axisymmetric expansion flow. For the compression flows, 
this model does not represent the integral length scale, because of the 
decrease in the dissipation in compression flows. 
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F. Testing Constants in Modeled Dissipation Equation 


The following is a method of deriving the model constants in the 
dissipation equation from the data. For homogeneous turbulent flows, 
the model equation for the dissipation rate is Eq. (7-2). In the iso- 
tropic compression flow, Eq. (7-2) can be rewritten: 

■111) - r 2c i + 3c 3 ) if) + c 2 (7 ‘ 10) 

e 

While in incompressible axisymmetrlc expansion flow, it can be rear- 
ranged to 


k d e 

2 W 


“ c , 17J + c 5 


(7-11) 


e 

We can use the simulation results to evaluate the validity of the model. 
The turbulence kinetic energy, dissipation, and production rate are 
obtained from the full- turbulence simulations. The time derivative of 
the dissipation is obtained by spline fitting the dissipation and dif- 
ferentiating the result; the derivatives obtained in this way contain 
considerable uncertainty. Despite the uncertainty, we should be able to 
discern trends; however, quantitative results need to be accepted with 
care. 


Equations (7-10) and (7—1 1 ) suggest that plots of - 
against P/e should be straight lines if model constants ar# indeed 

pure constants. Figures 7-21 and 7-22 show this kind of plots for iso- 
tropic compression and incompressible axisymmetrlc expansion flows. The 
high strain-rate simulations (runs SQF and EXQ) are not used here, due 
to unacceptable uncertainty caused by too few data points. As can be 
seen from Fig. 7-21, for isotropic compression, the slopes (-2Cj + 3 C 3 ) 
are the same for the various runs, but the intercepts (C 2 ) are not. 
We conclude that the model needs modification. Figure 7-22 shows that 
the situation is even worse for incompressible axisymmetrlc strain flow. 


Thus, the parameters in modeled dissipation equation cannot be pure 
constants if these data are to be fit. Dependence of the parameters 
on P/e is precluded by the linearity of Fig. 7-21. We considered the 
possibility that the model constants depend on Reynolds number. 
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However, the Reynolds number Is large enough that the model constants 
should not depend on this parameter. Furthermore, when the "constants" 
were plotted vs. Reynolds number, the variation of was opposite 
what one would expect. The remaining possibility is that the parameters 
are functions of the total strain ratio. 

To test this possibility, we re-sorted the data from the various 
runs according to the value of the total strain. For each total strain, 
a least squares fit to the data was made. Good fits were obtained, 
indicating that the model constants can be considered functions of total 
strain. k-e models that incorporate this dependence would require 
solution of several additional differential equations which describe the 
evolution of model constants* This procedure is cumbersome and not good 
modeling practice. 

We have now tested the k-e models as thoroughly as we can. The 
summarized results will be presented In the next section. 

G. Summary of k~e Model Testing 

We found that k-e models perform well in low strain-rate flows in 
which the production of energy is of about the same order of magnitude 
as the dissipation, i.e., the turbulence is in an "equilibrium" state. 
When the strain is so strong that the flow structure is out of equilib- 
rium, k-e models do not perform well. Furthermore, the model turbu- 
lence length scale does not represent the Integral length scale well in 
the compression cases. Finally, the "constants” in the modeled dissipa- 
tion equation cannot be pure constants. 

Urns k-e models do not perform well for high strain rate. A tur- 
bulence model capable of predicting flows accurately with a wide range 
of strain rates is needed. In the next section we shall propose a 
three-equation turbulence model which meets this need. 
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7.3 A One-Point-Closure Three-Equation Turbulence Model 

A. Motivation 

In two-equation (k-e) models, the turbulence kinetic energy de- 
termines the turbulence velocity scale, while the dissipation plays two 
roles: it is the rate of destruction of turbulence kinetic energy, (Eq, 
(7-1)), and it determines the turbulence length scale (Eq. (7-7)). The 
dissipation cannot do both jobs simultaneously in high strain-rate flows 
in which the flow structure is out of equilibrium. To decouple the dis- 
sipation and the length and time scales and introduce the minimum of 
additional complexity, a one-point-closure, three-equation turbulence 
model is proposed. In particular, a model equation for a turbulence 
time scale (t) is introduced. It is to be solved together with the 
dynamic equations for the turbulence kinetic energy (k) and its dissi- 
pation rate (e). At present, this model equation is limited to homo- 
geneous flows. The method of evaluating the constants of this model is 
discussed and tests of the performance of the model are presented below. 

B. Homogeneous Isotropic Decay Flow 

He begin by looking at the model as applied to Isotropic decay. 
The three-equation model consists of three differential equations. For 
isotropic-decay, the exact dynamic equation of the turbulence kinetic 
energy (Eq. (7-1)) reduces to: 


dk 

■at " e 


(7-12) 


The model equation for the dissipation rate is modified to include the 
new time scale: 


de _ _ e 

dt t 


(7-13) 


where t is the new turbulence time scale; this equation replaces Eq. 
(7-2). 

We need a model equation for t to close this set of equations; 
for isotropic decay it can be derived as follows. 
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The turbulence kinetic energy decays according to the power law 
(Reynolds, 1976) : 

k ~ (t - t o j“ n (7-14) 

The exponent (n) is reasonably well established as 1.2 at high Reynolds 
number; t Q is an effective origin. The dissipation rate (e) then 
decays as 


e 


dk 

dt 


- fc <J 


-n-1 


Then Eq. (7-13) requires 


T 


or, in differential form: 


nTT i* “ V 


dr _ n _ _5 
dt " n + 1 “ 11 


(7-15) 


(7-16) 


(7-17) 


Equations (7-12), (7-13), and (7-17) give the correct behavior of the 
turbulence in isotropic-decay. Of course, the added complication is 
unnecessary in this case. Also the constant in Eq. (7-17) should be a 
function of Reynolds number; this function could be evaluated using the 
data of Shirani et al. (1981). 


C. Return to Equilibrium 

When strain is applied to the flow, all the turbulence quantities 
are modified. In particular, the turbulence time scale (r) is pushed 
away from equilibrium, (cf. Eq. (7-16)). After the strain is removed, 
the turbulence tends to return to an equilibrium state. Hie simplest 
modification of the time-scale equation that will accomplish this (and 
be dimensionally correct) is : 


dr 

dt 




(7-18) 


Die form of the second term is selected because, in isotropic turbu- 
lence, ex/k ■ 6/11. Defining z « ex/k and manipulating Eqs. (7-12), 
(7-13), and (7-18), we obtain 
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If ■ |l* + SJ t'-TTJ <7 - 19) 

In isotropic-decay, Eq. (7-19) becomes trivial. Thus, in the (z,dz/dt) 
phase plane (6/11,0) is an equilibrium point. 

If z is perturbed from its equilibrium point in the direction of 
larger z and 1 + C 5 >0, then dz/dt will be positive and z will 
be driven away from the equilibrium state; if z is displaced in the 
opposite direction and 1 + C5 >0, z will also be driven away from 
the equilibrium state. Therefore, 

C 5 < - 1 (7-20) 

is required to assure return to equilibrium. 


D. Isotropic Compression Flow 


When turbulence is subjected to strain, the large eddies interact 
with the mean flow and extract kinetic energy from the mean motion. For 


isotropic-compression flows, the mean-deformation tensor U 




is 


shown in Eq, (2-26). The turbulence energy production is P 

2 — DIL 

k U .In the presence of homogeneous isotropic compression, the 

3 k,k 

three-equation model can be written 


dk 

dt 


DIL 


(7-21) 


P e 

d e e _ DIL e 

dt t 4 k 


(7-22) 


dr 

dt 


-1+ C .!!- 
11 5 1 k 


_6 

11 


I 4* G S x 

' 6, ISO ISO 


(7-23) 


The additional terms reflect the effects of the mean flow on turbulence 
quantities. Note that one could use 1/x in place of e/k 1° the last 
term of Eq. (7-22). The choice used in Eq. (7-22) was tried first and 
worked; other authors may wish to try the other possibility. 
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The three model constants (C^, C 5 , and C 6>IS0 ) can be evaluated 
from the results of lsotroplc-compresslon simulations. A method of 
accomplishing this follows. The turbulence kinetic energy, dissipation 
rate, and production used here are obtained from full-turbulence simula- 
tions. The time derivative of the dissipation is obtained by spline 
fitting the dissipation and differentiating the result. The slope of the 
curves in Fig. 7-21 suggests « 1. If this value is accepted, the 
turbulence time scale (f) is the only remaining unknown in Eq. (7-22) 
and can be evaluated. 

The model constants €5 and Cg i$q* are then evaluated from Eq. 
(7-23). The turbulence kinetic energy, dissipation rate, and mean 
strain rate (Sj§q) are obtained from full-turbulence simulations. The 
turbulence time scale (t) is obtained from Eq. (7-22) as described 
above, and its time derivative is obtained by spline fitting the time 
scale and differentiating the result. The model constants C 5 and 
Cg j S q are the remaining unknowns in Eq. (7-23) and can be evaluated. 
Keeping in mind that model constant C 5 must be less than -1 (Eq. 
(7-20)), we find that C 5 - - 1.1 and Cg jIS0 - - 0.5. 

The performance of the three-equation model for isotropic compres- 
sion with four different compression rates (runs SQF through SQI) is 
shown in Figures 7-23 through 7-26, respectively. The predictions of 
the k-e models are plotted for comparison. The evolution of the model 
time scale for each simulation is shown in Figure 7-27. As can be seen, 
the three-equation model is better than both the LSW and Reynolds' mod- 
els at all strain-rates. This is not surprising, because we fit the 
model constants to this flow. 
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where S is a positive constant. Ihe turbulence energy production 

through mean motion is P,. » - R D J , . The three-equation model in 

INC ij i,j 

the presence of homogeneous incompressible mean strain can be written 


dk 

dt 


P INC ~ e 


(7-25) 


d c e A . P INC e 

dt T 1 k 


(7-26) 


dr 

dt 


_5 

11 


C i— - -—) 
5 1 k 11 j 


+ C St 
6 ,AXI T 


(7-27) 


Ihe constants multiplying the strain terms are allowed to differ from 
those used in the isotropic compression case. The model constant C 5 
was determined to be -1.1 in the preceding section. The remaining two 
unknown model constants (Cj and Cg ^ axj) can be evaluated from the 
results of homogeneous, incompressible, axisymmetric , expansion flow 
simulations . 

The method used to evaluate the model constants Cj and Cg is 
similar to that presented in the preceding section. The turbulence 
kinetic energy, dissipation rate, and production are obtained from full- 
turbulence simulations. The time derivative of the dissipation Is ob- 
tained by spline fitting the dissipation and differentiating the result. 
This leaves two unknowns (Cj and x) to be determined in Eq. (7-26). 
It was found that, for incompressible axisymmetric expansion flows, the 
"production of dissipation" predicted by k-e models is too weak and 
the dissipation is underpredicted in high strain-rate flows. This means 
that Cj must be larger than 1.44 (the value in the LSW model) for this 
type of flow. By numerical experiments, we found Cj - 2. Equation 
(7-26) then provides x» Given the evolution of turbulence kinetic 
energy, dissipation rate, strain rate (from full-turbulence simula- 
tions), turbulence time scale (from Eq. (7-26)), and the time derivative 
of the turbulence time scale (from numerical differentiation), model 
constant Cg of Eq. (7-27) is found to be -2. 
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Since the model constants are evaluated from the "actual produc- 
tion", one needs an accurate Reynolds stress model to have a satisfac- 
tory prediction. The Boussinesq Reynolds stress model does not work 
well when the strain rate is high; its drawbacks were discussed in the 
preceding section. Therefore, a Reynolds stress model has been devel- 
oped for this flow. It will also be applied to the one-dimensional 
compression flow in the next section. 


F. Reynolds Stress Model 

An exact transport equation for the components of the Reynolds 
stress tensor can be obtained by multiplying Eq. (2-10) by uj , adding 
the result to the equation obtained by switching the subscripts 1 and 
j , and averaging the result : 


»R 


— + U R 


at 


k ij,k 


P + T - D + J 

ij ij ij 


(7-28) 


where pR - p u'u' is the Reynolds stress tensor. 

P ij ij 

The first term on the right side of Eq. (7-28), called the "produc- 
tion tensor", is the creation of Reynolds stress from the mean flow. 


here S 

1 TT 


^j 

1 7T 


- (R..s 


+ R ik S kiJ + 


Ik kj jk ki 


+ R jk n kiJ 


(7-29) 


— (U + U j is the mean strain rate tensor and Ujj 

Ij 2 i , j j , 1 J 


— I U. - U. ■) is the mean rotation 
2 1 i,j j,i J 


tensor , 


l ij 


is the "transfer tensor", 





u j.i J 


(7-30) 


It has zero trace in an incompressible turbulence field and is suppos- 
edly responsible for intercomponent redistribution of the Reynolds 
stress. It has the form of a correlation between the fluctuating pres- 
sure and the fluctuating strain rate. 


Djj is the "dissipation tensor". 



2v u' 


i,k j,k 


(7-31) 


It dissipates turbulent kinetic energy through viscous action. 


103 



The tensor is the diffusive flux of R^ , 


1 


J ** — i p'u' 6 + p*u' fi I + u'u'u’ - vR (7-32) 

ijk - i °jk p j °ik J ijk v ij,k V ' 

P 

By definition, homogeneous flows have no convection or diffusion, 
and Eq. (7-28) reduces to 


dR 


it 


dt 


P ij + T ij ' °ij 


(7-33) 


Pjj is explicitly determined from the Reynolds stress, but models are 
needed for T ij and . 

We first address the pressure-strain terms. An exact equation for 
the fluctuating pressure, p', can be derived by taking the divergence 
of Eq. (2-10). We get 


1 


P 







(7-34) 


The source term in this Poisson . equation can be split into two parts, 
giving rise to two components of the pressure: a "rapid" part p^ 
given by 


and a "slow" part 




given by 




(7-35) 


1 


P 





(7-36) 


Note that the "rapid" pressure, p^\ involves the mean deformation 
explicitly. Hence the imposition of a mean deformation Immediately 
changes the "rapid” pressure. Solving Eq. (7-35) for p^ 1 ^ (see 
Reynolds, 1984), we get 


where 



2U A 4 

p.j ijpq 


(7-37) 
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(7-38) 


M ijpq 



dk 


Models for **ijpq have been proposed by various modelers (Hanjalic and 
Launder, 1972; Launder, Reece, and Rodi, 1975; Hanjalic and Launder, 
1976; Lumley, 1978; and Lumley, 1979). The particular model we use was 
proposed by Reynolds (1984) and is similar to the others. It has the 
following form for divergence-free turbulence field. 


M 


ijpq 


t! Vp, * 15 ' V Jq + SV • V s y b p. 


+ VVJ, 


+ 6. b. + « b. + 6 b 

iq jp jp iq jq ip 


J + lJ 


1 4 


3 V 6 Pq b ij 


(7-39) 


where ■ 2k, b^ , defined in Eq. (6-2), is the Reynolds stress 
anisotropy tensor, 6 is the Kronecker delta, and is the model 
constant. The rapid part of the pressure-strain term is then modeled as 



2 2 

— S S 6 + 4A S b - 6A | b S + b S 

5 ij 15 kk°ij 1 kk ij l'- ik kj jk ki 


“ 3 b nm S nra 6 ij ” ^3 + 3 V ^ik^kj + hjk^ki- 1 ) 


(7-40) 


By insisting that this model agree with rapid distortion theory, the 
constant Aj is found to be -2/7 (Reynolds, 1984). For slower mean 
deformation rates, Reynolds (1984) added quadratic terms in bjj to 
M^pq. This complicates the model. 

In order to keep the model simple without losing the essence of 
modeling philosophy, we model T^j^ by using linear terms in b^j only 
(Eq. (7-40)), but allow the coefficient Aj to be a function of the 
ratio of production to dissipation to compensate. 


In one-dimensional compression or axisymmetric expansion flow, 
T^ is diagonal. Furthermore, since T^ * T^^ and T^ is 

traceless, the entire tensor can be expressed in terms of Tjj\ for 
which the model (Eq. (7-40)) becomes: 


(11 8 

T ii “ - 5 ks + 4 V 3R n " 2k J s 


(7-41) 
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where S is the mean strain rate of Eq. (7-24). Full- turbulence simu- 
lation (Lee and Reynolds, 1985) provides the actual rapid pressure- 
strain (TjpJ, turbulence kinetic energy (k), Reynolds stress (R^), 
and mean strain rate (S). These data allow us to evaluate the model 
coefficient Aj as a function of the ratio of energy production to dis- 
sipation (P/e “ (3Rn “ 2k)S/e), We found 


Aj - - 0.34 + 0.12 exp(-0.3 P/ e ) (7-42) 


Figure 7-28 shows how Eq. (7-42) fits the full-turbulence simulation 
data. As can be seen, Reynolds’ value (Aj ■ -2/7) is good for moder- 
ate and high strain rates (runs EXP and EXQ), to which rapid distortion 
theory may be applied. 


It is customary to model the combination of the slow pressure- 
strain term and dissipation anisotropy tensor together (Lumley and 
Newman, 1977; Lumley, 1978; Lumley, 1979; Rogallo, 1981; and Reynolds, 
1984). This tensor is defined as 


r ij 


■=■ p<2> ^i,3 + u i,i J ' * W 3 J 

0 


(7-43) 


where Dj^ - 2e. is modeled by 


♦ij * A o b ij (7 " 44> 

where bjj is the anisotropy tensor of Rjj and Aq is a model con- 
stant. In this flow, $21 is the only independent component of • 
We evaluated A Q from $21 and b ^ 2 obtained from full simulations. 
Figure 7-29 shows the evolution of A Q for flows with a wide range of 
strain rates. As can be seen, there is a considerable scatter. How- 
ever, $21 18 not important in moderate and high strain-rate cases 

(runs EXP and EXQ), because the mean strain rate is so high that mean 
flow dictates the evolution of turbulence. We chose A^ ■ 1 to fit the 
low strain-rate flow. Note that this choice might not give accurate 
modeling of the relaxation following removal of strain. For further 
results on modeling this term, see Lee and Reynolds (1985). 


Combining Eqs. (7-33), (7-29), and (7-43), the modeled Reynolds 
stress transport equation for homogeneous flow is written as 
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(7-45) 


dR 


ij 


dt 


n A J l > 
P, + T 

ij ij 


^ij ” 3 e6 ij 


where P^ , j ^ are given by Eqs. (7-29), (7-40), and (7-44), 
respectively. When this Reynolds stress model is applied to incompress- 
ible axisymmetric flow, it further reduces to the following form (bear 
in mind that only dRjj/dt is needed). 



_ A (1) 

P n T n " e *n 


2 

3 


e 


where 


P 


11 



S 


( 1 ) 

11 


- kS + 4A 1 ^3R n - 2k j S 


Aj = - 0.34 + 0.12 exp(-0.3 P/e) 
P « ^3R 1X - 2k j S 

♦ ll “ W 2k ~ l/3j 


A 


o 


1 


(7-46) 


Note that R22 * R33 ■ (2k - Rjj)/ 2 in this type of flow. 

The performance of this model (Eqs. (7-25), (7-26), (7-27), and 
(7-46)) for axisymmetric expansion flow at three strain-rates (runs EXQ 
through EXO) are shown in Figures 7-30 through 7-32, respectively. The 
predictions of the k-e models are plotted for comparison. The evolu- 
tion of the model tine scale for each simulation is shown in Figure 7- 
33. In high and moderate strain- rate flows (runs EXQ and EXP), the 
better performance of the k-g-t model can be attributed largely to the 
introduction of the new Reynolds stress model. In low strain-rate flow 
(run EXO), the k-g-r model is slightly better than k-g models, and 
there is no distinct difference between the LSW and Reynolds' models. 
Again, the k-g-t model is tuned to the full simulation results for 
this flow. It is expected to perform well. 


107 



G. One-Dimens lonal Compression Flow 


The mean-deformation tensor for one-dimensional compression is 



0 

0 

0 



(7-47) 


where S is negative. One-dimensional compression can be considered as 
a combination of the isotropic-compression and axisymmetrlc expansion. 
Formally, 



Both the dilatation and divergence-free parts of the mean flow contrib- 
ute to the energy production. The k~e~T model has the following form 


dk 

It 


P INC + P DIL * e 


(7-49) 


de _ _ e . r P INC e c P DIL g 

dt T C 1 k + C 4 k 


(7-50) 


dx 

dt 


11 


+ C_ 


1 k 


_6 

11 


^ + ^6, ISO 3 T + 


C — t 
6.AXI 3 


(7-51) 


where p xuc “ (2k - 3R 11 )S/3 and P DIIi - -2kS/3. The model constants 
obtained from the isotropic-compression and axlsymmetric expansion flows 
will be used without modification. They are given in Table 7.2. The 
Reynolds stress model tuned to Incompressible axisymmetric expansion 
flow (Eq. (7-46)) is also employed here. Since the dilatation part of 
the production, P DIL , does not contribute to the rapid pressure-strain 
tensor, the production in the expression for the model coefficient Aj, 
Eq. (7-42), should be ^XNC* This model (Eqs. (7-46), (7-49), (7-50), 
and (7-51)) reduces to those for the two preceding cases in the approp- 
riate limits. 
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Table 7.2 


The Values of the Model Constants In the k-g-t Model 


C 1 

2.0 


c 4 

1.0 


C 5 

-1.1 


C 6,IS0 

-0.5 

(for isotropic compression flow) 

C 6,AXI 

-2,0 

(for axisymmetric expansion flow) 


The three-equation model (Eqs. (7-49) through (7-51)) and the Rey- 
nolds stress model (Eq. (7-46)) were used to predict one-dimensional 
compression simulations with four strain rates (runs ODB through ODE). 
The predictions are shown in Figures 7-34 through 7-37. The predictions 
of the two-equation models are plotted for comparison. The evolution of 
the model time scale for each sllmulatlon is shown in Figure 7-38. In 
the rapid distortion case (run ODB), the k-g-t model is superior to 
the k-e models, due to the introduction of an accurate Reynolds stress 
model. In moderate strain-rate cases (runs ODC and ODD), the three- 
equation model reflects the behavior of the turbulence better than do 
the two-equation models. In slow strain-rate case (run ODE), there is 
little difference between the k-e-t and k-e models. 

It is encouraging that a model derived from other flows can be used 
to predict these flows and, in particular, that the effects of the 
various strains appear to be additive. 


H. Summary 

In this section we summarize the one-point-closure, three-equation, 
turbulence model. For homogeneous flows, this model is 


dr 

dt 


de 

dt 


dk 

dt 


P + p 
r INC r DIL 


_ £ • P INC e P DIL e 

x 1 k 4 k 


(7-49) 

(7-50) 


if + C s 


11 - 


+ C S T+C S T (7-52) 
6, ISO ISO 6,AXI AXI 
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where 


INC 


* - R 


1^1,: 


-An 


3 k,k ij 




DIL 


- — kU 
3 k,k 


S is the 


strain rate, and the are the model constants tabulated in Table 

7.2. Equation (7-49) is an exact equation for turbulence kinetic 
energy. 


Equation (7-50) is a model equation for the dissipation rate (e). 
The first term on the right side of Eq. (7-50) is the destruction term. 
The second and third terms are the "production of dissipation" due to 
the Incompressible and dilatation parts of the mean flow. The model 
constants Cj and are evaluated from the results of axisymmetric 
expansion and isotropic-compression flow simulations, respectively. It 
was found that - 2 and - 1 . 

Equation (7-52) describes the evolution of the turbulence time 
scale (t). The first term on the right gives the correct behavior of 
the turbulence quantities in isotropic-decay flow. The second term is a 
re turn-to-equi librium term, and the rest of the terms represent the in- 
fluence of the mean strain on the turbulence time scale. Stability 
analysis shows that model constant Cj must be less than -1. We chose 
C 5 to be -1.1 from numerical experiments. The new model constants are 
Cg jgQ “ -0.5 and Cg * -2. Three independent strain flows are 
needed to complete the model of the effect of strain on the turbulence 
time scale. In this work we covered two independent strains (isotropic 
compression and incompressible axisymmetric expansion), so one addi- 
tional term may be needed in Eq. (7-52) to produce a complete model. 
The effects of the strains appear to be additive. One more building- 
block flow such as axisymmetric contraction or plane strain is needed to 
calibrate the remaining model constant in the turbulence time-scale 
equation. 


A new Reynolds stress model was developed by fitting the simulation 
results for incompressible axisymmetric expansion flow. It has the form 


dR 


ij 


dt 


( 1 ) 2 
P ij + T ij - e *ij-3 e *ij 


(7-45) 
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(7-29) 



'ij ” " i R ik S kj + R jk S ki ^ + ^ R ik fi kj + R jAi J 


1 tr 


s ij “ + 

n ij ‘ 


q2 1 ! s u - if s i* s « + 4 a iVij - 6 vv« + v* 


3 b nm S nm 6 ij J " h + 3 V^ik^kj + b jk ft ki J 


(7-40) 


'ij 


- 0.34 + 0.12 exp(-0.3 P^/eJ 

R ij 6 lj 
3 


(7-42) 


♦ij 


A b b ij 


(7-44) 


Model constants ^ and Aj are tuned to fit Incompressible axisymmet- 
ric expansion flow; they were also used in predicting one-dimensional 
compression flow, and their performance was good. This model has not 
been tested in flows with rotation, and so should be used with caution 
in shear flows. 

The performance of the k-e~T model is slightly better than that 
of the k-e models in low strain-rate flows. In high strain-rate 
flows, the flow structure is out of equilibrium and the k-e~x model 
reflects the changing physics much better than do the k-e models. 
Some engineering flows occur in the range in which the differences of 
the models are significant. In summary, a one-point-closure, three- 
equation turbulence model which accurately calculates four types of 
flows — isotropic decay, isotropic compression, axisymmetric expansion, 
and one-dimensional compression flows- — for a wide range of strain rates 
has been developed. 
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TIME TIME 

Fig. 7-8(c). The prediction of R 

Fig. 7-8 (b) . The prediction of the dissipation models for run EXP. 

rate by k-e models for run EXP. 
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TOTAL STRAIN TIME 

Fig. 7 - 11 . The prediction of the turbulence pig. 7-12(3). The prediction of the turbulent 

length scale by Reynolds' model for kinetic energy by k-e models 

three axisymmetric expansional f or run qdb. 

strain simulations. 
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Fig. 7-14(b) . The prediction of the dissipation 7-14(c). The prediction of 

rate by k-e models for run ODD. models for run ODD. 
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Pig. 7-18(b). The evolution of a longitudinal Pig. 7-19(a). The evolution of q J /e In incora- 

lntegral length scale In Iso-* presslble, axlsymmetrlc expan- 
tropic compression simulations. slonal strain simulations. 
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Fig. 7-20 (b) . The evolution of the longitudinal Fig. 7-21. A plot of - (k/e 2 )(de/dt) against 

integral length scale In the P/e for Isotropic compression simu— 

compression direction for one- lations. 

dimensional compression Simula- 
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Pig. 7-26(b) • The prediction of the di68lpa— Pig. 7—27, The evolution of the modeled turbu— 
tion rate by k-e-r and k-e lence time scale In isotropic corn- 

models for run SQI. pression simulations. 




Fig. 7-28. A plot of the model coefficient Aj Fig. 7-29. The evolution of the model coeffi- 

as a function of P/e and its fit- cient A^ in Incompressible, axi- 

ting in incompressible, axisymmet- symmetric expansional strain flow, 

rlc expansional strain flow. 
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Fig. 7-30(a). The prediction of the turbulent Fig. 7-30(b). The prediction of the dissipation 

kinetic energy by k-e**T and rate by k-e-t and k-e models 

k-e models for run EXQ. ^ or run 
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Fig. 7-32(a). The prediction of the turbulent Fig. 7-32 (b) . The prediction of the dissipation 

kinetic energy by k-e-t and rate by k-e-t and k-e models 

k-e models for run EXO. for run EXO. 
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Fig. 7-35(a). The prediction of the turbulent 
Fig. 7-34 (c) . The prediction of Rjj by k-e-x kinetic energy by k-e-x and 

and k-e models for 1 run ODB. k-e models for run ODC. 
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Fig. 7-35 (b). The prediction of the dissipation 

rate by k-e-T and k-e models Fig. 7-35(c). The prediction of Rjj by k-e— 

for run ODC. and k-e models for run ODC. 




0 . 5 



147 


Fig. 7-36(a) . The prediction of the turbulent Fig. 7-36(b). The prediction of the dissipation 

kinetic energy by k-e—r and rate by k-e-t and k-e models 

k-e models for run ODD. for run ODD. 
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-38. The evolution of the modeled 
turbulence time scale In one- 
dlmenslonal compression flow. 



Chapter VIII 

CONCLUSIONS AND RECOMMENDATIONS 


Conclusions and recommendations distilled from the results and 
discussion presented in the preceding chapters are given in this chap- 
ter. 


8.1 Conclusions 

Numerical methods for solving the three-dimensional, time- 
dependent, Navier-Stok.es equations for homogeneous turbulence undergoing 
isotropic and one- dimensional compression have been presented along with 
results obtained from such simulations. The simulations were carried 
out at low Reynolds number; the limitation on this parameter derives 
from computer resource availability. The results obtained for fast 
compression rates have been compared with rapid distortion theory, and 
good agreement was obtained. 

The ratio of the turbulence time scale to the Imposed mean-flow 
time scale was found to be the most important single parameter in these 
flows. When this ratio is large, the flow is immediately affected by 
the mean strain in a manner similar to that predicted by rapid distor- 
tion theory; i.e,, the turbulence kinetic energy and its dissipation 
rate increase, while the integral length scales immediately decrease. 
When this ratio is small, the flow initially retains the character of 
decaying isotropic turbulence; only after the strain has been applied 
for a long period does the flow begin to reflect the effect of mean 
strain. In these flows, the Kolmogorov length scale, which is the size 
of the eddies in which viscous dissipation occurs, decreases rapidly 
with decreasing total strain, due to the density increase that accom- 
panies compression. This means that these flows have more turbulence 
energy in the small scales than would be expected on the basis of the 
models usually used to simulate these flows. 

Results from isotropic compression, incompressible axisymmetric 
expans ional strain, and one-dimensional compression simulations were 
used to test one-point-closure, two-equation (k-e) turbulence models. 
It was found that k-e models perform well at low strain rates. When 
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the strain is so strong that the flow structure is out of equilibrium, 
these models do not perform well. The parameters in the model dissipa- 
tion equation cannot be pure constants if the model is to fit the 
computer-generated data. It was also found that the model turbulence 
length scale, which is supposed to represent the integral length scale, 
does so very poorly in these flows. In compressed turbulence, the rap- 
idly decreasing Kolmogorov scale caused by the density increase drives 
the flow out of equilibrium. The dissipation can no longer satisfactor- 
ily fulfill the two functions that k-e models require of it; it cannot 
simultaneously be both the rate of destruction of turbulence kinetic 
energy and the length scale determiner. 

To decouple the dissipation and the time scale, a one-point, three- 
equation turbulence model including a model equation for the turbulence 
time scale is proposed. The time-scale equation contains terms which 
account for the decay of isotropic turbulence, the return to isotropic 
turbulence after a strain has been removed, and the effects of various 
types of strain. The other two equations of the model govern the evo- 
lution of the turbulence kinetic energy and its dissipation rate and are 
similar to those in k-e models. The essential change is that the time 
scale in the destruction term in the modeled dissipation equation is 
replaced by the new time scale. It is encouraging that the effects of 
various strains appear to be additive, i.e*, only few building-block 
flows are needed to calibrate the model constants. The new model 
accurately calculates four types of flows — isotropic decay, isotropic 
compression, axisymmetric expansion, and one-dimensional compression 
flows — for a wide range of strain rates. 

8.2 Recommendations 

As the speed and memory si2e of supercomputers advance, it will 
become possible to do simulations of turbulence on a 256 x 256 x 256 
grid. This will allow us to study the evolution of turbulence which has 
an inertial subrange in the spectrum and may eliminate the questions 
associated with the low Reynolds numbers of the present simulations. 

Not all of the flows needed to calibrate the new three-equation 
turbulence model were simulated in this work. The present code (with 
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slight modifications) could be used to simulate all the building-block 
flows needed to calibrate the new model. The model should also be 
tested by applying it to the simulation of inhomogeneous flows. 

Swirl is deliberately generated in the intake stroke by design 
inlet port shape and orientation and persists through to top-dead-center 
(TDC) of the compression stroke, at which point the mean flow pattern 
approximates solid-body rotation. The effect of swirl on turbulence is 
not well understood. It could be studied by modifying the present code 
to include a Coriolis force term in the Navier-Stokes equations. 

Finally, homogeneous turbulence undergoing expansion could be simu- 
lated with the present code with minimum modifications. Results from 
such simulations could help explain the turbulence behavior during the 
expansion stroke in an internal combustion engine. 
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Appendix 
TABULATED DATA 


Data tabulated in this appendix were obtained from computed turbu- 
lence fields described in Chapters V and VI, Fields from the following 
runs are included. 


Run ID 

Strain Type 

Mean-Deformation Tensor U 

SQF 

Isotropic compression 

Eq. (2-26) with L Q - 0.3, V p - -5.6 

SQG 

•# it 

" " " L o -1.0, V p «-1.0 

SQH 

** #• 

" " " L -0.3, V - -0.06 

o ’ p 

SQI 

it ii 

L q - 0.3, V p = -0.012 

EXQ 

Incomp, axisym. exp. 

Eq. (7-24) with S - 35.87 

EXP 

It II It 

• S - 3.587 

EXO 

•I •• •• 

" " " S - 0.3587 

ODB 

1— D compression 

Eq. (2-30) with L Q ■ 0.3, V p - -24.3 

ODC 

•t it 

L q - 0.3, V p - -1.29 

ODD 

ii «• 

" " " L q - 0.3, V p - -0.26 

ODE 

ii <i 

L q » 0.3, V p - -0.052 


The following data are presented for each field : 


t 

S 

TS 

v 


time 

mean strain rate 
total strain “ exp( 
kinematic viscosity 



S(t’)dt’J 


with units [T] 

ir 1 ] 

dimensionless 

[LV 1 ] 


Strain Type 

S 

TS 

V 


V 

IRH1H 

O 1 

Isotropic compression 

^(t) • ^ + y t 


v(t) - v(0)(TSr ,A 


0 p 

o 


Incompressible, axi- 

constant 

exp(S(t-t Q )) 

constant 

symmetric expansion 


t Q - 0.58789 



V 

L +V t 

r» 7 

1-D compression 

" L +V t 

. £ E. 

L 

v(t) - v(0)(TS) U,/ 


o p 

o 
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uju^ ^ no summat * on ) 


with units [L 2 T“ 2 ] 


D. . “ 2v ut .ul . (no summation on i) 


q z * Rjl + R22 + R33 “ twice the turbulent 

kinetic energy 


[i 2 f 2 j 


t - <D„ + D 22 + D33>/2 * dissipation rate 


U.V 3 ] 


r h 

L,. “I Q, . (x_) dx (Lv * half the computational box) with units [L] 
ij * m J n 13 m m n 


longitudinal 


lateral 


“ -(=^M 

you */axj)y 


integral length scales in x m ~direction 


(no summation) 


Taylor microscales 


The following data are presented for incompressible, axisymmetrlc expan- 
sion simulations only: 

dR;i/dt with units [L 2 T“ 3 ] 


[L 2 r 3 ] 


^ ~p a> »i,i 

p 


[l 2 t" 3 1 


(D n - 2e/3) 


[L 2 T -3 ] 
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